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ABSTRACT 


Based  on  potential  flow  theory,  a  formulation  is  given  for  three- 
dimensional  fully  cavitating  flow  with  a  Riabouchinsky  model.  The  model 
is  nonlinear  and  the  location  of  the  free  surface  of  the  cavity  is  not 
known  priori.  Therefore,  an  iterative  procedure  is  used  to  locate  the 
free  surface  boundary.  The  employment  of  a  trial-free-boundary  approach 
effectively  reduces  the  fully  nonlinear  model  to  a  linear  one,  and  the 
solution  at  each  iteration  is  obtained  by  means  of  the  finite  element 
method  (FEM) .  Examples  studied  were  fully  cavitating  flow  past  flat 
plates  in  a  water  tunnel.  Results  are  given  for  pure  drag  flows  past 
circular  and  elliptic  plates  and  a  lifting  flow  past  a  circular  plate. 

Theree-dimensional,  20-node,  quadratic  isoparametric  finite  elements 
are  used.  The  locations  of  the  mid-side  nodes  on  the  free-surface  just 
off  the  edge  of  the  plate  are  found  to  have  significant  effects  on  the 
separation  condition  which  must  be  enforced  while  keeping  the  Jacobian 
of  the  isoparametric  transformation  nonsingular.  The  free  surface  is 
shown  to  be  a  characteristic  surface  and  this  leads  to  the  development 
of  a  weighting  scheme  for  assigning  the  free  surface  potential  in  the 
iterative  scheme  the  shifting  of  the  free  surface.  A  new  three-dimensional 
"local"  iterpolation  scheme  for  the  location  of  the  mid-side  nodes  is 
also  developed  for  the  lifting  flow  configuration. 

The  various  algorithms  were  tested  against  known  analytic  solutions 
as  well  as  published  numerical  and  experimental  results  in  two-dimensional 
and  axisymmetric  flows.  Results  from  several  pure-drag  three-dimensional 
geometries  are  presented  and  the  cavity  shapes  and  the  tunnel  wall  effects 
are  examined.  Then  lifting  flow  results  are  given. 

Because  of  the  change  in  flow  boundary  conditions  at  the  separation 
edge  and  the  failure  of  the  FEM  to  resolve  these  conditions  accurately, 
the  ability  of  the  numerical  solution  to  maintain  a  constant  pressure  over 
the  entire  cavity  decreases  as  the  three  dimensionality  of  the  free  surface 
increases.  However,  the  present  procedure  produces  absolutely  stable 
iterations  and  shows  no  sign  of  drifting  of  the  free  surface.  It  is 
found  that  satisfaction  of  a  tangent  separation  condition  of  the  free 


surface  from  the  flat  plate  body  is  crucial  for  the  stability  of  the  iterative 
procedure.  Grid  refinement  in  both  the  streamwise  and  transverse  directions 
reduces  the  computational  error.  While  free  surface  movement  between 
iterations  is  a  useful  convergence  criterion,  a  flow-rate  balance  between 
upstream  and  downstream  cross-sections  appears  not  to  be  a  good  criterion. 
Finally,  alternate  formulations  which  may  reduce  the  difficulties 
encountered  are  discussed. 
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1.  INTRODUCTION 


1.1  OVERVIEW 

The  analytical  study  of  wake  and  cavity  flows  began  in  the  19th 
century  when  Kirchoff  introduced  an  idealized,  inviscid-flow,  infinite- 
length-cavity  model  with  free  streamlines  and  solved  the  accompanying 
problem  via  the  conformal-mapping  technique.  Since  then,  numerous 
researchers  have  investigated  different  aspects  of  the  problem,  both  ana¬ 
lytically  and  numerically.  The  development  of  underwater  missiles,  high 
speed  hydrofoils  and  propellers  accelerated  the  research  efforts  on  the 
analysis  of  and  design  for  cavity  flows.  Models  were  subsequently  formed 
to  investigate  flows  with  finite  cavities.  The  Riabouchinsky  model,  the 
re-entrant  jet  model,  the  open  wake  model  and  the  vortex  model  are  some  of 
those  being  used  (see,  e.g.,  Wu,  1972). 

Cavity  flow  problems,  as  represented  by  potential  flow  models,  are 
nonlinear.  The  problems  are  complicated  by  the  unknown  shape  and  location 
of  the  free  boundary.  For  two-dimensional  cavity  flows,  methods,  such  as 
the  hodograph  mapping  which  is  derived  from  the  theory  of  analytic  functions, 
can  be  employed  to  transform  the  problems  onto  complex  planes  where  the 
problem  is  solved  by  singular-integral-equation  techniques.  The  evaluation 
of  these  functional  equations  is  very  difficult  in  general  because  of  the 
nonlinearity.  But,  by  and  large,  the  two-dimensional  problem  for  steady 
state  flows  has  been  solved  (Tulin,  1963;  Street  and  Larock,  1968;  Wu,  1968). 

Though  efforts  have  been  made  to  solve  the  three-dimensional  case, 
there  is  still  no  known  exact  solution  or  generally  usable  and  fully  non¬ 
linear  model  for  the  analysis  of  the  three-dimensional  fully  cavitating  flow 
(see  the  review  of  numerical  methods  for  solution  of  three-dimensional  cavity 
flow  problems  by  Street,  1977).  He  points  out  that  the  techniques  for  the 
numerical  solution  of  three-dimensional,  fully  cavitating  flows,  based  on 
linearized  methods  and  the  method  of  matched  asymptotic  expansions,  are 
presently  available.  However,  he  also  points  out  that  these  techniques 
rely  on  the  use  of  two-dimensional  characteristics  of  the  flow  applied  strip- 
wise  in  the  three-dimensional  flow  field  and  that  the  area  of  fully  nonlinear 
methods  is  still  developing. 
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Among  Che  linearized  models  are  Widnall's  (1966)  lifting-surface 
model,  which  was  found  Co  compare  unsatisfactorily  with  experimental  data 
by  Tsen  and  Guilbaud  (1974) .  The  disagreement  was  credited  to  the  non¬ 
linearity  of  the  physical  problem.  Nishiyama  (1970),  with  the  assumption  of 
large  aspect  ratios,  was  able  to  reduce  the  lifting-surface  problem  to  a 
lifting- line  problem  but  achieved  similar  results  as  Widnall  (1966) .  Jiang 
and  Leehey  (1977)  combined  the  essence  of  these  two  theories  and  produced 
calculations  that  appear  to  be  the  best  within  the  limitations  of  the 
linearized  theory. 

Furuya  (1975) ,  using  two-dimensional  nonlinear  streamline  theory  in 
the  near  field  and  Prandtl’s  lifting  line  theory  for  the  far  field,  and 
Leehey  and  Stellinger  (1975),  using  higher  order  perturbation  theory,  have 
gone  beyond  the  linearized  models  and  developed,  not  linear,  but  yet  not 
fully  nonlinear  models.  Neither  model  accounts  for  wall  effects,  although 
Furuya's  (1975)  model  can  handle  the  motion  of  foils  beneath  free  surfaces 
in  water  of  infinite  depth.  A  three-dimensional  fully  nonlinear  model  is 
needed  for  the  study  of  hydrofoils  with  finite-aspect-ratio  and  the  inter¬ 
pretation  of  the  results  from  water  tunnel  tests.  It  was,  therefore,  the 
goal  of  this  research  to  develop  a  computational  procedure  for  the  simula¬ 
tion  of  three-dimensional  flow  about  fully  cavitating  hydrofoils. 

Numerical  techniques  which  are  commonly  used  in  free  surface  flow 
computations  are  reviewed  in  Section  1.2.  Section  2.1  describes  the  physi¬ 
cal  formulation  while  Section  2.2  gives  the  finite  element  (FE)  approximation 
of  the  three-dimensional  fully  cavitating  flow  problem.  The  computation  of 
flow  rates  and  areas  in  the  FE  representation  is  described  in  Section  2.3. 
Section  2.4  describes  the  characteristic  of  the  free  surface  which  leads  to 
the  development  of  a  weighting  scheme  described  in  Section  3  together  with 
the  iterative  procedure  for  the  location  of  the  free  surface.  Section  4 
then  gives  the  implementations  and  tests  of  the  algorithms  while  computa¬ 
tion  results  are  given  in  Section  5.  Section  6  contains  a  discussion 
errors  and  alternative  approaches.  Section  7  combines  the  conclusions  and 
suggestions  of  possible  future  work. 
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1.2  REVIEW  OF  NUMERICAL  TECHNIQUES 

There  are  several  major  categories  of  numerical  techniques  which  can 
be  applied  to  free  surface  flows.  There  are  finite-difference  and  finite- 
element  methods.  Amongst  the  finite-difference  methods  one  finds  formula¬ 
tions  in  which  the  problem  is  solved  in  a  stream  function  and  velocity 
potential  space  or  in  which  the  solution  is  found  from  the  distribution  of 
singularities  on  the  free  surface  of  the  flow  (normally  called  a  boundary 
integral  technique) .  The  application  of  all  these  techniques  to  both  two- 
and  three-dimensional  free  surface  flow  problems  is  described  below. 

1.2.1  Two-Dimensional  Methods 

Mogel  and  Street  (1974)  used  a  finite-difference  technique  and  solved 
a  problem  of  flow  past  a  disk  in  a  water  tunnel  by  employing  a  Riabouchinsky 
model  and  imaging  to  achieve  an  exact  formulation.  Irregular  finite- 
difference  stars  were  employed  along  the  curved  free  surface  boundary.  A 
very  refined  finite-difference  grid  was  also  used  in  the  neighborhood  of  the 
separation  point  on  the  flat  plate.  The  solution  of  the  finite-difference 
equations  for  the  velocity  component  was  achieved  by  successive  over¬ 
relaxation.  Because  the  velocity  components  were  employed  as  the  unknowns, 
the  solution  required  simultaneous  satisfaction  of  the  Laplace  equation  for 
each  velocity  component  and  of  the  nonlinear  free  surface  boundary  condi¬ 
tions.  The  solution  converged  and  gave  reasonable  results.  However,  the 
cost  for  the  two-dimensional  case  was  subsequently  found  to  be  equal  to  the 
cost  for  a  solution  of  similar  accuracy  for  a  three-dimensional  disk-in- 
water-tunnel  flow  by  the  finite-element  method  (see  Sec.  5). 

Brennen  (1969)  combined  a  finite-difference  technique  with  a  mapping 
to  stream  function — velocity  potential  space  to  solve  the  axisymmetric  flow 
past  a  circular  disk  in  a  circular  water  tunnel  using  a  Riabouchinsky 
model.  A  difficulty  with  the  extension  of  this  technique  to  more  general 
bodies  is  that  the  solution  is  Inverse  because  the  stream  function — velocity 
potential  space  solution  must  be  mapped  back  to  physical  space  to  obtain 
the  final  geometry.  Brennen 's  results  are  used  for  comparisons  with  the 
finite  element  solution  obtained  from  the  three-dimensional  model  later  in 
the  present  work. 
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White  and  Kline  (1975)  developed  a  general  method  for  the  solution  of 
turbulent,  separated,  axisymmetric  flows.  In  their  general  development 
they  employed  a  new  boundary  integral  technique  for  the  potential  flowfield 
outside  the  boundary  layers  developing  in  their  diffuser  flows.  Unlike 
grid  techniques,  the  boundary  integral  method  requires  computation  of  un¬ 
knowns  only  over  the  flow  boundary.  The  solution  technique  obtains  the 
potential  flow  from  the  numerical  solution  of  Green's  third  identity  and 
iteratively  solves  for  the  correct  location  of  an  initially  unknown  free 
surface.  In  principal  this  technique  can  be  extended  to  general  three- 
dimensional  flows.  Interestingly,  a  very  similar  technique  was  developed 
in  1953  by  Armstrong  and  Dunham  (1953) .  They  applied  a  vortex  sheet  singu¬ 
larity  to  free  streamline  flows,  thereby  generating  an  integral  equation 
whose  solution  formed  part  of  an  iterative  procedure  to  locate  the  initially 
unknown  free  surface.  They  solved  several  axisymmetric  cavity  flows  in  an 
infinite  fluid  field  and  obtained  reasonable  agreement  with  experiment. 

The  application  of  the  finite  element  method  to  free  streamline  flows 
was  first  made  by  Chan,  et  al.  (1973  a,b) .  They  showed  that  the  finite 
element  method  could  be  used  to  solve  for  the  velocity  potential  in  two- 
dimensional  and  axisymmetric  ideal  fluid  flows  involving  a  free  surface. 
Sarpkaya  and  Hiriart  (1975)  applied  a  similar  technique  to  solve  the  free 
streamline  problem  of  flow  in  curved  jet  deflectors. 

1.2.2  Three-Dimensional  Methods 

Jeppson  (1972)  developed  an  inverse  formulation  for  three-dimensional 
flows  in  which  he  defined  a  velocity  potential  and  two  additional  functions 
which  defined  the  flow  paths.  By  changing  the  conventional  roles  played  by 
the  variables  of  the  problem,  he  converted  a  free  surface  with  an  unknown 
position  in  physical  space  into  a  plane  of  known  position  in  the  inverse 
stream  function — velocity  potential  space.  The  technique  is  very  similar 
to  that  of  Brennen  (1969)  sited  above  and  suffers  the  same  disadvantage, 
namely  that  the  shape  of  curved  solid  bodies  cannot  be  prescribed  in  advance 
in  physical  space. 

Larock  and  Taylor  (1976)  applied  finite-element  techniques  to  solve 
the  jetflow  from  a  circular  pipe  and  orifice  under  the  influence  of  gravity. 


The  resulting  flow  is  slightly  three-dimensional,  namely,  the  jet  remains 
essentially  circular  but  droops  under  the  action  of  gravity  creating  a 
three-dimensional  flowfield.  They  employed  the  velocity  potential  as  the 
dependent  variable  and  quadratic  isoparametric  finite  elements  in  the  flow- 
field  discretization.  Again,  as  in  other  free  streamline  problems,  an 
iterative  procedure  had  to  be  used  to  determine  the  free  surface  location 
whose  precise  location  is  initially  unknown.  They  found  "convergent"  solu¬ 
tions  when  the  jet  is  not  too  far  from  being  horizontal  (i.e.,  at  high 
Froude  number).  However,  they  reported  that  "a  number  of  small  shifts  (of 
the  free  surface  between  two  consecutive  iterations)  can  accumulate  to  a 
nonnegligible  effect."  Plausible  reasons  for  the  inability  of  their  pro¬ 
cedure  to  converge  at  low  Froude  numbers  and  for  the  drift  of  their  solu¬ 
tion  are  advanced  later  herein  as  the  present  work  is  described. 

In  summary,  the  works  of  Mogel  and  Street  (1974),  Brennen  (1969), 
and  Larock  and  Taylor  (1976)  are  of  special  interest.  Although  developed 
only  for  two-dimensional  flow  and  based  on  a  finite  difference  scheme,  Mogel 
and  Street  (1974)  provided  the  basic  framework  for  the  present  study  and 
development  of  a  three-dimensional  finite  element  model.  Brennen  (1969) 
found  good  agreement  between  his  solution  and  experiments.  Hence,  his 
problem  was  solved  herein  by  the  finite  element  method  as  a  check  of  the 
present  three-dimensional  model.  Larock  and  Taylor  (1976)  provide  an  illus¬ 
tration  of  the  application  of  the  finite  element  method  to  a  three- 
dimensional  problem  and  gave  us  insight  to  problems  to  be  expected  and  key 
features  to  be  incorporated  into  a  model. 
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2.  FORMULATION  FOR  FULLY  CAVITATING  FLOWS 


2.1  PHYSICAL  MODEL  DESCRIPTIONS 

One  common  consideration  which  enters  into  the  development  of  cavity 
flow  models  is  the  problem  of  the  closure  of  the  cavity.  In  the 
Riabouchinsky  model,  the  closure  of  the  cavity  is  achieved  by  a  mirror 
image  of  the  foil  in  the  downstream  direction  as  shown  in  Figure  2.1, 
with  section  C  being  the  plane  of  symmetry.  In  the  present  study,  the 
Riabouchinsky  model  centered  within  a  water  tunnel  of  uniform  geometric 
cross-section  is  used.  It  is  assumed  that  both  ends  of  the  water  tunnel  can 
be  extended  far  away  from  the  foils  so  that  a  uniform  horizontal  velocity 
can  be  maintained  at  the  entrance.  Because  of  the  symmetry  in  the  flow 
geometry,  only  one-eighth  (for  the  pure  drag  case)  or  one-quarter  (for  the 
lifting  case)  of  the  flow  field  need  to  be  considered. 

The  flow  is  assumed  to  be  incompressible,  steady  and  irrotational. 
Hence,  it  is  governed  by  the  velocity  potential  <j>  which  satisfies  the 
Laplace  equation,  in  Cartesian  coordinates, 

Ijt  + 

3x2 

within  the  flow  field  and 

u 

V 

w 

where  u,  v  and  w  are  the  components  of  the  velocity  in  the  x,  y,  and  z 

222  1/2 

directions,  respectively.  Their  magnitude  q*(u+v+w)  .  If  the 

effect  of  gravity  is  neglected,  the  Bernoulli  equation  becomes 


+  IJL  »  0 
2  2 
3y  3z 


3$ 

3x 


3(1) 

3y 


3d) 

3z 


(2.1) 


(2.2) 


(2.3) 


1  2 

P  +  p  q  ”  constant 

where  P  is  the  pressure  and  p  is  the  fluid  density.  Since  the  free 
surface  is  also  a  constant  pressure  surface,  from  the  Bernoulli  equation. 
Equation  2.3, 


u+v+w  =  q  »  constant 

on  the  free  surface  and  the  change  of  velocity  potential  <j> 
surface  along  a  streamline  is  given  as 


(2.4) 

on  the  free 


3<fr 

3s 


q 


(2.5) 


Also  the  free  surface  can  be  described  by  the  streamline  conditions,  namely, 


^  =  4*  (2.6) 

U  V  w 

which  mean  that  the  velocity  vector  is  tangent  to  the  free  surface. 

Because  of  the  symmetry  of  the  flow  model,  the  plane  of  symmetry 
(section  C,  Figure  2.1)  of  the  Riabouchinsky  model  is  also  an  equipotential 
surface.  Without  loss  of  generality,  the  potential  there  can  be  set  to 
zero.  With  a  prescribed  uniform  horizontal  upstream  flow  (with  velocity 
)  and  no  flux  conditions  at  solid  boundaries  (tunnel  walls  and  the  foil) , 
the  problem  is  well  defined.  However,  the  location  of  the  free  surface  is 
not  known  a  priori.  An  iterative  procedure  is  described  in  Section  3  for 
determining  the  location  of,  as  well  as  the  constant  pressure  condition  on, 
the  free  surface.  As  a  result  of  using  such  a  trial-free-boundary  method, 
only  a  linear  system  of  equations  has  to  be  solved  at  each  iteration. 

The  most  important  nondimens ional  parameter  for  cavity  flows  is  the 
cavitation  number  a  defined  as 


a 


(2.7) 
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where  the  subscript  c  denotes  the  cavity.  Equality  of  cavitation  numbers 
is  an  indication  of  the  dynamic  similarity  of  two  geometrically  similar 
cavity  flows.  Using  Equation  2.3,  one  can  show  that 


2 

c_ 

2 

00 


-  1 


(2.8) 


The  force  on  the  foil  is  given  by 


/ 


A 


(P  -  P  )  dA 
P  c 


(2.9) 


where  subscript  p  denotes  the  foil  (a  plate  herein)  and  A  is  the  area 
of  the  foil.  From  Equations  2.3  and  2.9, 


*f  f' s  ^ 

A 


(2.10) 


For  the  pure  drag  problem,  F  is  the  drag.  For  a  plate  at  an  angle  of 
attack  a  the  drag  D  is  given  by 


D  ■  F  sin  a 


and  the  lift  L  is  given  by 


L  *  F  cos  a  . 


(2.12) 


The  drag  and  lift  coefficients,  and  respectively,  are  given  as 

follows : 


i  pui s 


(2.13) 


and 


\  3 


(2.14) 


where  S  is  defined  as  the  largest  possible  projected  area  of  the  foil. 

It  is  often  convenient  to  normalize  and  on  the  velocity  q^  on 

the  free  surface.  Then, 

C?  -  CD/(l+a)  -  :  °  2  (2.13a) 

2  S 

and 

C*  -  C  /(1-HJ)  -  .  L ,  (2.14a) 

2  S 

where  the  reslationship  between  a,  q^  and  is  given  by  Equation  2.8. 

2.2  FINITE  ELEMENT  MODEL 

The  finite  element  method  (FEM)  was  originally  developed  in  the  1950's 
in  the  field  of  structural  engineering  and  used  because  of  its  efficiency  in 
stress  analysis  of  large  structural  systems.  Since  then,  the  FEM  has  become 
one  of  the  commonly  used  methods  for  the  approximate  numerical  solution  of  a 
wide  spectrum  of  engineering  problems.  One  of  the  advantages  of  the  finite 
element  approximation  is  its  simplicity  and  efficiency  for  handling  irregu¬ 
lar  boundaries  and  variation  in  mesh  sizes  which  are  necessary  to  the  solu¬ 
tion  of  free  surface  problems. 

In  the  FEM,  the  solution  domain  is  discretized  into  subregions,  called 
elements,  which  are  defined  by  their  node  points  (often  located  at  the 
corners  and  mid-sides  of  elements).  The  FEM  then  represents  the  problem  solu¬ 
tion  as  a  continuous  function  by  using  approximation  or  interpolation  func¬ 
tions  within  »ach  element.  These  functions  are  defined  in  terms  of  the 
solution  values  at  the  finite  set  of  node  points.  The  solution  varies  locally 
within  each  element  according  to  the  form  determined  by  the  interpolation 
function;  these  normally  are  prescribed  via  a  summation  of  a  set  of  shape  (or 
basis)  functions.  Partial  differential  equations  are  then  transformed  to 
integral  equation  formulations  which  are  subsequently  replaced  by  a  well- 
posed  set  of  algebraic  equations  for  the  values  at  the  node  points. 
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2.2.1  Finite  Element  Approximation 

As  described  in  Section  2.1,  the  fully  cavitating  flow  problem  is 
nonlinear  because  of  the  unknown  location  of  the  free  surface.  However, 
the  nonlinear  solution  can  be  achieved  by  solving  a  series  of  linear  prob¬ 
lems  through  the  use  of  an  iterative  procedure  which  is  described  in  Sec¬ 
tion  3.  It  is  the  linear  problem  associated  with  each  iteration  which  is 
solved  by  means  of  the  FEM.  In  particular  the  velocity  potential  is 
determined  for  a  given  configuration  then  the  result  is  used  to  determine 
a  more  accurate  free  surface  shape.  As  the  application  of  the  FEM  to  gen¬ 
eral  field  problems,  i.e.,  the  determination  of  <J>  ,  is  well  documented 

(see,  e.g.,  Huebner,  1975),  only  a  brief  description  is  given  here. 

The  solution  domain  &  is  discretized  into  the  subregion  elements. 

Each  element  is  defined  by  node  points  as  shown  in  Figure  2.2.  The  contin¬ 
uous  unknown  function  <j>  is  approximated  by  the  interpolation  function  deter¬ 
mined  by  the  unknown  values  (f^  at  all  the  node  points  in  the  domain.  The 
local  variation  of  the  solution  within  each  element  is  determined  by  the 
shape  functions  because  the  solution  is  expressed  as  a  linear  combination  of 
the  products  of  the  values  of  the  shape  functions  and  the  ,  defined  at 
each  node  point  of  the  element,  i.e., 


<pe  -  E  N.  <p.e  (2.15) 

i-1  X  1 

where  the  superscript  e  denotes  values  for  element  e  .  The  subscript  i 

til 

denotes  values  defined  for  the  i  node  point,  n  is  the  total  number 
of  node  points  for  each  element  and  are  the  shape  functions;  <J>^  is 

a  subset  of  <p ^  . 

Shape  functions  are  classified  by  two  characteristics,  viz.,  first, 
the  degree  of  variation  within  each  element,  such  as  linear,  quadratic, 
etc.,  and,  second,  the  order  of  continuity  of  the  solution  function  <|> 
between  elements,  e.g.,  C°  continuity  to  ensure  the  continuity  of  the 
solution  function  between  boundaries  of  elements ,  continuity  to  ensure 

the  continuity  of  the  solution  function  as  well  as  its  first  derivatives 
between  boundaries  of  elements,  and  so  on.  However,  all  shape  functions 
satisfy  the  following  condition: 
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• 


(2.16) 


The  total  number  of  node  points  n  per  element  is  fixed  by  the  choice  of 
shape  functions.  The  choice  of  the  shape  function  for  the  present  study 
is  addressed  later  in  this  section. 

The  linear  problem  to  be  solved  at  each  iteration  is  the  Laplace 
equation  for  the  velocity  potential  ,  namely, 


k(%)  ’  0 


in  Q, 


(2.17) 


(summation  notation  is  used)  with  boundary  conditions 


(a)  <p  -  $(xi)  on  Sx 

(b)  l.  +  g(x.)  =0  on  S, 

i  i  < 


(2.18) 


(2.19) 


where 


are  the  direction  cosines  of  the  outward  normal  to  the  surface 


(Figure  2.2).  This  problem  is  actually  solved  by  application  of  the 
variational  (Ritz)  approach,  which  uses  the  variational  form  of  the  differ¬ 
ential  equation  rather  than  the  differential  equation  itself.  Thus,  the 
solution  <p  minimizes  the  functional 


m%)' 


dx^  + 


(2.20) 


It  can  be  shown  that  there  exist  a  set  of  values  ,  defined  at  all  the 
node  points  of  a  finite  element  mesh,  which  yield,  through  the  basic  func¬ 
tions,  the  best  approximation  to  the  continuous  solution  <j>  of  the  differ¬ 
ential  equation.  It  can  also  be  shown  that  this  best  approximation  is 
unique  (Prenter  1975,  p.  196).  If  the  total  functional  is  taken  to  be  the 
sum  over  all  elements,  minimization  of  the  functional,  Equation  2.20,  leads 
to  the  following: 
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i 

\  3z  )\ 

(2.21) 


The  validity  of  the  summation  on  the  righthand  side  of  Equation  2.21  is 
assured  by  the  right  choice  of  shape  functions.  Zienkiewicz  (1971,  Section 
3.2)  gives  two  criteria  to  be  satisfied  by  the  chosen  shape  functions.  The 
shape  functions  chosen  herein  meet  these  criteria;  in  particular,  as  only  the 
function  $  and  its  first  derivatives  appear  in  (2.21),  a  C°  element  and 
its  shape  functions  are  needed. 


(2.22) 


From  Equations 

2.15  and 

2.21,  one 

obtains,  for  each  element, 

„e 

e  e 

0 

Kij 

3  * 

the  elemental  " 

'stiffness 

"  matrix 

Kij 

is  given  by 

i re  s 

Kij 

rn_h 

3Nj 

3N  3Nj  \ 

J  \  3x 

_  p 

3x  3y 

3y 

(2.23) 


and  the  forcing  term  due  to  the  flux  boundary  condition  on  is  given  by 


>i  •  /  «»l 


dS„ 


(2.24) 


Note  that,  if  none  of  the  node  points  of  an  element  is  on  S£  ,  the  forcing 

e 

term  vanishes.  The  global  equation  is  obtained  by  summing  over  the 

elemental  equations,  i.e., 
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(2.25) 


Equation  2.25  is  a  set  of  linear  algebraic  equations  from  which  the  <|> 
can  be  obtained.  The  global  stiffness  matric  is  positive  definite  and 

symmetric.  Both  direct  and  iterative  methods  for  solving  this  large  sparce 
linear  system  are  reviewed  in  Appendix  2. 

2.2.2  Isoparametric  Transformation  and  Numerical  Integration 

In  order  to  handle  irregular  boundaries  and  variation  in  mesh  sizes, 
a  special  class  of  finite  elements,  the  isoparametric  finite  elements,  is 
used.  The  idea  of  using  elements  with  curved  sides  seems  to  have  originated 
with  Taig  (1961)  and  was  generalized  by  Irons  (1966)  and  Ergatoudis,  et  al. 
(1968) .  The  isoparametric  transformation  uses  a  set  of  shape  functions  to 
map  an  irregularly  shaped  element ,  which  is  defined  in  physical  space  by 
means  of  its  node  points ,  onto  a  regular  cube  defined  in  a  local  coordinate 
system  as  shown  in  Figure  2.3. 

Herein,  a  set  of  quadratic  shape  functions  is  employed.  Accordingly 
in  three-dimensions,  each  isoparametric,  C°  ,  finite  element  has  20  node 
points.  The  transformation  is  given  explicitly  as 


x  = 


20 

Z  N  .x, 
i=l  1  1 


y  * 


20 

Z  N  .y . 

•  i  i  i 
i=l 


(2.26) 


z  = 


20 

Z  N.z, 
,  ,  i  i 
i=l 


where  IT  *  N^(5,n,i;)  are  the  shape  functions.  For  isoparametric  elements. 


employed  in  the  coordinate  transformation  (Equation 


the  shape  functions 

2.26)  are  the  same  as  those  employed  to  represent  the  unknown  field  variables 
(see  Equation  2.15).  The  quadratic  shape  functions  appropriate  for  the 
20-node  "serendipity"  element  (Huebner,  1975)  used  here  are: 
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(2.27) 


for  corner  nodes,  i.e.,  5^  -  =  4^  ■  +  1  : 

»i  -  jd  +  «,)(!  +  nnt)  (l  +  44^(5^  +  nni  +  44*  -  2) 


for  midside  nodes  at  "  0  ,  n  *  5  -  +  1 

N±  -  ^(1  -  C2)(i  +  nn.Ki  +  44*) 


for  midside  nodes  at  ni  *  0,  *  +  1 

Ni  -  i(i  +  cc±)  (1  -  n2)d  +  44*) 


for  midside  nodes  at  c.  =  0,  £.  ■  n,  *  +  1 

i  i  i  — 

n*  =  J(i  +  55±)  (l  +  nn*)  (i  -  42) 


These  satisfy  the  constraint.  Equation  2.16. 


The  essence  of  the  transformation  is  contained  in  the  Jacobian  matrix 
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functions  are  related  as 
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It  can  also  be  shown  that 


dxdydz  -  (j|  d£dnd£ 


(2.30) 


However,  this  transformation  is  valid  if  and  only  if  the  Jacobian  |j|  is 
nonsingular  in  an  element.  If  the  transformation  is  consistent,  then  the 
sign  of  the  Jacobian  jjj  is  the  same  in  all  elements.  It  is  well  known 
(see,  e.g.,  Strang  et  al.  1973,  3.3)  that  all  interior  angles  at  the 
corners  of  an  element  should  be  bounded  between  0  and  tt  ,  and  the  midside 
nodes  must  lie  in  the  middle  half  of  the  sides  of  the  element  to  guarantee 
nonsingular  Jacobians.  It  is  also  true  that  unique  relations  between  the 
physical  (x,  y,  z)  and  local  (£,  r\,  £  )  coordinate  systems  do  not  exist 
for  elements  which  fold  back  upon  themselves,  i.e.,  those  with  singular 
Jacobians. 

Whenever  an  unique  transformation  exists,  the  integration  of  any 
function,  say  u(x,  y,  z) ,  defined  in  physical  space  can  be  performed  in 
isoparametric  space  (local  coordinate  system)  such  that 


x,y,z)  dxdydz 


+L  ±1  ±1 

ILL 


u(x(5,n,0 ,  y(5.n,s) , 


z(C,n,0)  * 


•  |J|  d^dndc 


(2.31) 


The  triple  integral  on  the  righthand  side  of  Equation  2.31  can  be  approxi¬ 
mated  by  means  of  a  Guass-Legendre  quadrature  formula  which  will  integrate 
exactly  a  polynomial  of  degree  2n  -  1  or  less  as  follows 


f*  2--1  d*  ’  *  “iP2n-l<:ti) 

J-l  i=l 


(2.32) 


where  x^  are  the  zeroes  of  the  Legendre  polynomial  of  degree  n  and 
are  the  weights.  For  a  triple  integral. 


J-l  J-l  J-l  n  n  n 

fff  u(5 »n,£)d£dr)dC  -  I  E  I  u^.n. ,?.  )  (2.33) 

J—  1  J—  1  J—  1  i=l  j»i  k=l  ^  ^ 
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The  locations  and  the  corresponding  weights  can  be  found  in  subroutine  Gauss 
described  in  the  Program  Manual  (Ko  and  Street,  1979). 

2.2.3  Implementation  of  Boundary  Conditions 

With  an  assumed  location  of  the  free  surface,  the  geometry  of  the 
flow  field  is  defined  and  can  be  discretized.  As  noted  in  Section  2.1, 
only  one-eighth  (for  the  pure  drag  case)  or  one-quarter  (for  the  lifting 
case)  of  the  total  flow  needs  to  be  considered  because  of  the  symmetry  of 
the  flow  geometry.  A  schematic  of  the  flow  field  employed  in  the  FEM  solu¬ 
tion  (for  a  pure  drag  case)  is  shown  in  Figure  2.4  for  the  purpose  of  illus¬ 
trating  the  boundary  conditions.  It  is  obvious  that  no-flux  conditions  are 
required  at  the  tunnel  wall,  the  plane  of  symmetry  and  the  surface  of  the 
foil  as  indicated  in  Figure  2.4.  Actually  the  free  surafce  is  also  a  no-flux 
boundary.  However,  the  no  flux  conditions  is  satisfied  only  if  the  free 
surface  is  at  its  true  location,  where  the  constant  pressure  conditions  is 
also  satisfied.  Therefore,  the  no-flux  condition  cannot  be  imposed  at  the 
free  surface  during  the  iterative  procedure,  but  the  final  solution  must 
satisfy  such  a  condition. 

% 

2. 2. 3.1  Flux  Boundary — Uniform  Flow  Upstream 
A  uniform  flow  U  is  imposed  at  some  reasonable  distance 

OO 

upstream  from  the  foil  or  plate.  Because  the  FEM  expresses  the  continuous 
solution  of  the  problem  in  terms  of  the  values  at  node  points,  flow  distrib¬ 
uted  over  the  surface  of  an  element  has  to  be  represented  by  the  equivalent 
nodal  values  specified  at  the  node  points  which  define  the  surface. 
Zienkiewlcz  (1971,  p.  168)  has  given  some  distribution  values  for  the  allo¬ 
cation  of  a  unit  uniform  surface  load  for  rectangular  surfaces  when  the  mid¬ 
side  nodes  are  placed  at  exactly  the  1/2  (quadratic  shape  function)  or  1/3 
(cubic  shape  function)  points  along  the  sides  of  the  elements.  However,  for 
surfaces  with  curved  sides  or  where  the  mid-side  nodes  are  not  at  symmetri¬ 
cal  locations,  one  has  to  evaluate  an  integral  for  the  equivalent  nodal 
values  E^  due  to  a  uniform  surface  load  as  follows : 

Ei  «  ff  UooN1  dxdy  -  dxdy  -  UJ^  (2.34) 
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where  is  the  equivalent  nodal  value  at  the  i  node,  Um  is  the 

uniform  surface  load  (in  the  present  case,  being  the  uniform  flow  upstream), 
is  the  shape  function  and  dx^y  is  the  weight  of  the  it*1 

node.  Also,  can  be  interpreted  as  the  area  of  influence  of  the  ith 

node  and  has  the  following  property: 


Z  u>j 


(2.35) 


where  n  is  the  total  number  of  node  points  on  the  surface  and  Ag  is  the 
area  of  the  surface.  The  evaluation  of  can  be  done  in  isoparametric 

space  as 


"i  "  //•  si 


n)  ( J |  d?dn  (2.36) 


where  ,  the  two-dimens ional  quadratic  shape  functions  as  used  in  the 
present  study,  are  given  as  (cf..  Equation  2.27  for  the  three-dimensional 
shape  functions) : 


for  corner  nodes:  i.e.,  ^  B  ±  “  +  1 

-  j(l  +  5?±)  (l  +  nn^  (CC±  +  nni  -  l) 

for  mid-side  nodes  at  =*  0  ,  =  +  1  : 

»i  -  |(1  -  S2)  (1  +  Tin±) 


(2.37) 


for  mid-side  nodes  at  *  0  *  a  +  1  : 


\  -  fa  +  cq)d  -  n2) 


and  the  Jacobian  of  the  transform  is  given  as: 


3x 

3x 

h. 

an 

3n 

(2.38) 


The  evaluation  of  the  integral  (Equation  2.36)  for  can  be  accomplished 

by  means  of  a  Guass-Legendre  quadrature  formula  similar  to  the  one  given 
by  Equation  2.33. 


2. 2. 3. 2  Boundary  Conditions  on  the  Free  Surface 


From  Equations  2.5  and  2.6,  it  follows  that  two  boundary  condi¬ 
tions  must  be  imposed  on  the  free  surface.  If  the  potential  at  the  equi- 
potential  plane  of  symmetry  (Section  AA,  Figure  2.4)  is  set  to  zero,  the 
potentials  at  the  node  points  along  any  assumed  streamline,  which  issues 
from  the  edge  of  the  plate  or  foil,  can  be  estimated  according  to  the 
constant  pressure  condition  (Equation  2.5)  by  integration  back  from 
Section  AA,  such  that  (referring  to  Figure  2.5), 


j 

i+1 


dS 


(2.39) 


where  the  subscript  i  refers  to  values  at  node  points  along  the  same 
streamline  and  the  superscript  j  refers  to  different  streamlines  on  the 
free  surface.  The  potentials  at  mid-nodes,  such  as  S^+'*'^  ,  are  inter¬ 
polated  from  the  values  at  the  adjacent  comer  nodes,  and  S^+^  .  To 

start  the  iterative  procedure,  q  is  estimated  by  use  of  the  upstream 
uniform  flux  and  a  mass  conservation  balance,  subsequently,  q  is  adjusted 
to  satisfy  the  separation  tangency  condition.* 

The  second  boundary  condition,  e.g.,  the  tangency  of  the  free  sur¬ 
face  to  the  velocity  vector  (Equation  2.6)  is  used  during  the  iterative 
procedure  to  create  new  (better)  free  surfaces  by  integration  downstream 
from  the  separation  surface  on  the  foil. 

2. 2. 3. 3  Effect  of  Mi-side  Node  Locations  at  Separation 

The  singular  properties  of  isoparametric  elements  have  been 
observed  and  studied  by  several  investigators  (e.g.,  Henshell  et  al.  1975, 


*In-flow  equals  out  flow  almost  exactly  once  the  solution  has  been  iterated 
a  few  times  and  a  mass  balance  is  no  longer  a  meaningful  test.  See 
Sec.  5. 
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Barsoum,  1976)  in  the  field  of  structural  analysis.  They  have  demonstrated 

-1/2 

the  possibilities  of  producing  strain  singularities  of  order  x  by 
placing  the  mid-side  nodes  of  two-dimensional,  quadratic,  isoparametric 
elements  at  the  quarter  point  along  an  element  side  from  the  corner  where 
the  singularity  occurs.  However,  it  can  be  shown  that,  as  a  consequence, 
the  Jacobian  of  the  isoparametric  transformation  is  singular;  therefore,  it 
may  generate  undesirable  behavior  near  the  singularity.  In  fact,  Hibbitt 
(1977)  showed  that  in  the  case  of  the  quadrilateral  element  used  by  Barsoum 
(1976) ,  the  strain  energy  of  that  element  is  unbounded  even  though  good 
numerical  results  were  reported  by  Barsoum  (1976) .  A  similar  procedure  is 
used  in  the  present  study  to  improve  the  behavior  of  the  solution  at  separa¬ 
tion;  however,  the  Jacobian  of  the  isoparametric  transformation  is  kept 
nonsingular  here. 

Consider  a  surface  of  an  element  just  off  the  edge  of  the  foil  as 
shown  in  Figure  2.6.  It  is  desirable  for  the  present  physical  realization 
to  have  (i)  the  edge  a  asymptotically  perpendicular  to  the  foil  at  the 
separation  point  S  ,  and  (ii)  the  edge  8  asymptotically  tangent  to  the 
foil  at  the  separation  (the  Kutta  condition).  From  the  isoparametric  map¬ 
ping  of  the  edge  a  ,  it  can  be  shown  that 


21  ~  4z9  +  3z13 
xi  “  4x9  +  3x13 


(2.40) 


To  satisfy  condition  (i) ,  one  needs  to  have  (dz/dx)|^  =  0  in  Equation 
2.40  which  leads  to 


'l  -  N  +  3z13 


-  21 


(2.41) 


It  is  clear  from  Equation  2.41  that  the  edge  a  has  zero  slope  at  the 
separation  point  S  if  the  mid-side  node  # 9  is  placed  at  the  quarter  point 
location  from  corner  node  #13  in  the  z-coordinate.  The  x-coordinate  of 
mid-side  node  #9  is  still  arbitrary.  However,  to  avoid  having  a  singular 


«.  .--tokaS* 


Jacobian  for  the  transformation,  the  x-coordinate  of  the  node  it 9  should  be 
kept  away  from  the  quarter  point  location.  Condition  (i)  also  ensures  the 
continuity  of  the  x-derivative  (x-component  of  the  velocity  vector  in  the 
present  study)  of  the  field  variable  as  one  moves  down  the  body  and  past 
the  separation  point  S  (across  the  element  boundaries)  and  onto  the  free 
surface.  This  character  is  essential  to  an  accurate  satisfaction  of  the 
Kutta  condition  (u  =  0)  on  the  horizontal  component  of  velocity  at  separa¬ 
tion. 

Similarly,  from  the  mapping  of  edge  (3  , 


dz  ® 
d*  13 


~3*13  +  *214  -  z15 
-3*13  +  4*14  -  *15 


(2.42) 


condition  (ii)  requires  Equation  2.42  to  be  infinite  at  the  separation  point 
S  .  This  is  achieved  by  setting  the  denominator  of  Equation  2.42  to  zero; 
hence , 

“3x13  +  4x14  '  x15  "  0 

/  X13  ~  x15  \  ,,  ... 

'  '  x14  "  X13  "  \  4  /  (2.43) 

Therefore,  as  far  as  the  local  shape  of  the  free  surface  is  concerned,  the 
tangent  separation  condition  is  satisfied  at  every  iteration  if  the  mid-side 
node  #14  is  placed  at  a  quarter  point  location  in  the  x-coordinate.  Similar 
to  the  case  for  condition  (i) ,  the  z-coordinate  of  the  mid-side  node  #14  can 
be  arbitrary  except  that  it  cannot  be  at  the  quarter  point  location  in  z 
because  then  the  Jacobian  |J  |  of  the  transformation  will  become  singular. 
However,  it  is  shown  later  that  the  effect  of  satisfying  the  "local"  separa¬ 
tion  condition  in  this  manner  is  not  significant.  One  can  conclude  that  the 
local  behavior  of  the  free  surface  within  the  element  at  separation  is  not  of 
paramount  importance;  however,  the  summed  behavior  (i.e.,  the  average  value 
of  the  velocity  at  a  node)  does  seem  to  be  crucial  as  shown  below. 


20 


'**./*  I  ••*‘1  '  f—  ■■wy*-***. 


2.3  COMPUTATION  OF  AREAS  AND  FLOWS 

In  the  course  of  the  iterative  procedure,  it  is  necessary  to  compute 
areas  of  arbitrary  plane  shapes,  flows  through  arbitrary  plane  cross- 
sections,  and  the  sum  of  squares  of  the  velocities  over  plane  foils  of 
arbitrary  shape.  The  area  of  the  foil  and  the  sum  of  squares  of  the 
velocities  over  the  foil  are  needed  in  the  computations  of  the  drag  and 
the  drag  coefficient.  The  trial  velocity  on  the  free  surface  at  Section  AA 
(Figure  2.4)  is  initially  estimated  from  the  total  flows  through  the  up¬ 
stream  and  downstream  cross-sections.  All  these  can  be  evaluated  in  a  simi¬ 
lar  fashion  in  the  isoparametric  space. 

For  areas  of  cross-sections  and  foils,  one  can  use  Equations  2.35  and 
2.36.  However,  the  area  of  any  surface  of  an  element  also  can  be  obtained 


+1  +1 


■  ffd*dy  -  f  f 


J  d£dn 


(2.44) 


-1  -1 


where  |j|  is  the  Jacobian  of  the  transformation  given  by  Equation  2.38. 
For  a  two-dimensional  isoparametric  element  with  n  node  points,  the 
Jacobian  |j|  of  the  transformation  can  be  computed  by 


J11J22  "  J12J21 


where 


n  3N 

Jl^lT  * 


n  3N 

<  si  “aT  ’ 


(2.45) 


n  3N, 


i  3n  ’ 


n  3N 

.1  yi 


and  Ni  is  the  shape  function  of  the  transformation. 


has 


Similarly,  for  a  flow  Q  through  arbitrary  cross-sections,  one 


which  is  needed  by  Equation  2.10  for  estimation  of  the  drag  and  lift 
coefficients.  Equations  2.44,  2.48,  and  2.49  can  be  evaluated  by  means  of 
a  Gauss-Legendre  quadrature  formula  in  two-dimensions  as  described  in 
Section  2.2.3. 
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2.4  CHARACTERISTICS  OF  THE  FREE  SURFACE 

As  mentioned  in  Section  2.1  an  iterative  procedure  is  needed  for 
determining  the  location  of,  as  well  as  the  constant  pressure  condition  on, 
the  free  surface.  The  constant  pressure  condition  is  satisfied  on  each 
streamline  at  each  step  (this  is  required  by  the  use  of  Equation  2.39); 
however,  the  constant  pressure  condition  does  not  need  to  be  satisfied 
across  the  streamlines  (i.e.,  the  pressure  on  each  streamline  may  be  dif¬ 
ferent  from  that  on  its  neighbor)  until  the  solution  has  converaged.  This 
characteristic  is  important  because  as  described  later,  it  suggests  a  stable 
iterative  procedure. 

The  location  of  the  free  surface  is  assumed  at  the  start  of  the  iter¬ 
ative  procedure.  Thereafter,  the  new  location  of  the  free  surface  is 
approximated  by  moving  the  present  free  surface,  based  on  the  streamline 
conditions  (Equation  2.6)  along  rays  of  constant  j  (Figure  2.5),  such 
that  the  new  free  surface  is  tangential  to  the  velocity  vectors  computed 
on  the  old  free  surface.  Velocities  on  the  free  surface  at  any  iteration 
are  known  functions  of  the  spatial  coordinates  (x,y,z  in  Cartesian  coordi¬ 
nate  system)  having  been  determined  from  differentiation  of  the  velocity 
potential  field  which  was  determined  for  a  given  geometry  prior  to  shifting 
of  the  free  surface. 

Generally,  any  given  surface  can  be  described  by 

F(x,y,z)  -  0  (2.50) 


or  specifically  by 


F  -  X(y,z)  -  x  -  0  (2.51) 

From  differentiating  Equation  2.51,  one  obtains 

dx  -  X  dy  +  X  dz  (2.52) 

y  z 

where  the  subscript  denotes  partial  differentiation  with  respect  to  that 
variable.  However,  for  the  free  surface  which  is  tangent  to  the  velocity 
field,  one  can  write  (cf.,  Equation  2.6) 
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dx 

dt 

iz 

dt 

dx 

dt 
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From  Equations  2.52  and  2.53,  one  finds 


vX  +  wX  - 
y  z 


u 


(2.54) 


as  the  equation  of  the  free  surface.  It  can  be  shown  that  (Street  1973, 
Chapter  9)  Equation  2.54  is  a  first-order,  quasi-linear  partial  differ¬ 
ential  equation  which  has  characteristics  given  by 


dx  =  d£ 
u  v 


dz 

w 


(.2.55) 


and  a  unique  solution  determined  by  the  passage  of  the  characteristics 
through  a  single  noncharacteristic  curve  on  the  surface.  Equation  2.55 
is  exactly  the  same  as  the  streamline  conditions  given  by  Equation  2.6. 
Therefore,  the  free  surface  is  a  characteristic  surface.  Furthermore,  the 
non-characteristic  curve  is  precisely  the  separation  curve  which  represents 
the  edge  of  the  foil.  Thus,  the  solution  to  the  problem  hinges,  first,  on 
an  accurate  integration  of  Equations  2.55  or  2.6  and,  second,  on  an  accurate 
representation  of  the  velocity  field,  in  particular,  the  conditions  at  the 
separation  curve.  As  we  demonstrate  below,  the  first  condition  is  easy  to 
satisfy;  the  second  is  not.  If  errors  creep  in  at  separation  the  exact 
characteristic  through  each  node  on  the  separation  curve  uniquely  repre¬ 
sents  different  flow  conditions.  Accordingly,  one  can  expect  progressive 
divergence  of  the  iterative  solution  (cf.,  Larock  and  Taylor,  1976,  where 
drift  of  the  solution  is  noted) . 

Only  when  the  solution  has  converged,  can  one  say  that  the  computed 
characteristic  surface  must  feel  the  same  pressure  everywhere  (i.e., 
constant  pressure  and  only  then  need  be  the  same  for  all  j  in 
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Figure  2.5).  [Note  that  the  magnitude  q  of  the  velocity  is  related  to 
the  pressure  by  Equation  2.3.1  Accordingly,  during  the  course  of  the 
iterative  procedure,  it  is  possible  to  adjust  the  q^  from  streamline  to 
streamline,  such  that  the  separation  velocity  at  each  node  tends  toward 
an  average  value  typical  for  all  nodes  on  the  separation  curve.  This  leads 
to  the  development  of  a  weighting  scheme  described  in  Section  3  and  tends 
to  keep  the  shifting  streamlines  on  the  same  characteristic  surface. 
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3.  ITERATIVE  PROCEDURE 

The  problems  which  remain  Co  be  solved  after  the  formulation  in 
Section  2  are  two.  First,  because  the  trial-free-boundary  method  is  used, 
the  problem  for  the  velocity  potential  <j>  is  reduced  to  a  well-posed 
linear  problem.  Under  the  finite  element  method  used  here,  a  linear  system 
of  algebraic  equations  is  all  that  must  be  solved  (See  Appendix  2) . 

Second,  however,  the  correct  location  for  the  cavity  free  surface 
must  be  obtained.  This  can  be  accomplished  by  an  iteration  procedure  in 
which  sequentially  one  prescribes  a  best  estimate  of  the  cavity  surface 
and  then  finds  the  corresponding  4>  values.  From  these,  velocities  are 
determined  and  they  give  the  basis  for  prescribing  a  new  free  surface 
location. 

3.1  OVERVIEW  OF  THE  ITERATIVE  PROCEDURE 

The  iterative  procedure  is  needed  to  determine  the  correct  location 
of  a  cavity  free  surface  which  is  tangent  to  the  flow  velocity  and  on  which 
the  constant  pressure  condition  is  accurately  satisfied.  The  step  by  step 
iterative  procedure  Is  described  in  this  sub-section.  Detailed  descriptions 
of  specific  schemes  are  given  in  the  following  two  sub-sections. 

The  procedure  consists  of  seven  steps  as  follows: 

Step  1:  Geometry  Setup 

With  given  shapes  and  sizes  of  the  plate  and  water  tunnel  and  an 
assumed  initial  location  of  the  free  surface,  the  solution  domain  is  well 
defined.  It  is  then  discretized  into  finite  elements  of  variable  sizes. 
Smaller  elements  are  used  in  regions  where  high  resolution  is  required,  for 
example,  near  separation  or  zones  of  rapid  changes  in  flow  direction. 

From  a  computational  point  of  view,  the  actual  discretization  of  the 
solution  domain  is  done  separately.  (The  mesh  generation  procedure  is 
described  in  Appendix  1.)  Then  the  coordinates  of  all  corner  nodes  together 
with  the  corresponding  nodal  numbers  become  inputs  to  the  actual  flow- 
computation  model.  The  mid-nodes  of  all  elements  are  interpolated  auto¬ 
matically  by  the  computer  program. 
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Step  2:  Boundary  Conditions 

An  uniform  horizontal  velocity  is  imposed  upon  the  upstream 

cross-section  of  the  water  tunnel  (cf.t  Section  2.2).  The  velocity 
potential  4>  on  the  free  surface  is  estimated  from  the  assumed  location 
as  described  in  Section  2.2.  At  the  solid  boundaries  and  the  planes  of 
symmetry,  the  no  flux  condition  is  imposed.  In  the  FEM,  this  natural 
boundary  condition  is  assumed  by  the  program  when  neither  <J>  nor  a  flux 
are  prescribed. 

Step  3:  Generation  of  Stiffness  Matrices 

With  the  established  coordinates  of  the  node  points  from  Step  1, 

0 

the  elemental  stiffness  matrices,  of  Equation  2.23,  are  generated. 

Because  of  the  great  amount  of  data  being  generated,  these  elemental  stiff¬ 
ness  matrices  are  usually  stored  on  magnetic  disks  or  tapes. 

Step  4:  Solution  for  the  Velocity  Potential 

The  velocity  potential  is  determined  by  solving  Equation  2.25.  How¬ 
ever,  the  global  stiffness  matrix  is  not  formed  explicitly.  Methods 

for  solving  this  large  sparse  linear  system  are  reviewed  in  Appendix  2. 

Step  5:  Computation  of  Velocity  Field 

The  fluid  velocities  are  calculated  by  differentiating  the  velocity 
potential  <J>  as 


3<P 

*7 


n 

l 

i=*l 


3N 

"5x7 


j  -  1,  2,  3 


(3.1) 


where  the  subscript  i  denotes  values  defined  for  the  ith  node  point, 
n  is  the  total  number  of  node  points  for  each  element  and  are  the 

shape  functions.  Equation  3.1  gives  the  fluid  velocities  in  three  direc¬ 
tions  at  all  node  points  of  element  e.  For  node  points  that  are  common 
to  more  than  one  element,  the  fluid  velocities  are  taken  to  be  the  alge¬ 
braic  average  of  the  contributions  from  all  adjacent  elements. 

Step  6:  Relocation  of  Free  Surface  and  Test  of  Convergence 

This  step  is  the  substantial  focus  of  this  work  and  many  details 
are  given  in  subsequent  sections.  To  avoid  the  evolution  of  a  divergent 
solution,  the  free  surface  location  is  actually  changed  only  when  the 
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magnitude  of  the  velocity  component  normal  to  the  plate  at  separation  is 
less  than  a  prescribed  value.  (This  leads  to  a  system  which  always  satis¬ 
fies  the  Kutta  or  tangent  flow  separation  condition.)  If  the  separation 
velocities  satisfy  this  condition  as  described  in  Section  3.3,  the  free 
surface  is  relocated  such  that  the  new  free  surface  is  tangent  to  the 
velocity  field  which  existed  at  corresponding  node  points  on  the  previous 
free  surface.  The  details  of  the  relocation  procedure  are  considered  in 
Section  4.1.  Next,  the  criterion  for  convergence  is  checked  (cf..  Section 
4.2).  If  it  is  satisfied,  the  iteration  is  terminated  and  the  solution  has 
been  obtained.  Otherwise,  the  mesh  in  the  near  vicinity  of  the  free  surface 
is  adjusted  and  the  elemental  stiffness  matrices  for  those  elements  which 
were  moved  are  replaced.  Since  only  those  elements  which  are  near  the  free 
surface  are  moved,  a  relatively  small  number  of  elemental  stiffness  matrices 
need  to  be  replaced. 

Step  7:  Boundary  Conditions  for  the  Free  Surface 

From  the  new  location  of  the  free  surface,  the  velocity  potential  <p 
on  the  free  surface  can  be  estimated  according  to  Equation  2.39.  However, 
the  separation  edge  represents  a  boundary  singularity  (Crank  and  Furzeland, 
1978)  at  which  two  types  of  boundary  conditions  meet,  namely,  the  no  flux 
condition  and  the  prescribed  velocity  potential  <(>  .  The  error  in  the 
solution  <}>  at  this  point  was  shown  by  Crank  and  Furzeland  (1978)  to  be  a 
function  of  nodal  placement;  a  similar  relation  can  be  expected  for  the  FEM, 
i.e.,  the  error  decreases  as  element  size  decreases  (see  Section  6). 

A  significant  and  supplementary  effect  may  enter  in  the  present  case, 
however.  In  the  true  three-dimensional  case  the  curvature  of  streamlines 
varies  along  the  separation  line  and  the  streamwise  rate  of  change  of  the 
potential  <j>  will,  therefore,  be  different  at  different  nodes  along  the 
separation  line.  Accordingly,  we  believe  that  the  boundary-condition- 
discontinuity  error  is  not  constant  along  the  separation  surface.  To  com¬ 
pensate  for  this  we  have  constructed  a  special  scheme  for  adjusting  the 
potential  <j>  on  individual  streamlines  (also  see  Section  6) . 

After  the  first  iteration  when  the  potential  <J>  is  calculated  from 
the  downstream  end  of  the  cavity,  we  alter  the  separation  values  of  the 
potential  and  then  prorate  that  value  to  other  nodes  along  the  surface  such 
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that  a  constant  velocity  is  obtained  on  each  streamline.  However,  this 
produces  a  different  constant  velocity  on  each  stream.1  ine.  As  the  solu¬ 
tion  converges  (in  the  sense  that  the  cavity  surface  does  not  move  signifi¬ 
cantly)  ,  one  hopes  that  the  final  velocity  differences  are  small.  As  is 
shown  in  Section  3.3,  where  the  weighting  scheme  is  described  in  detail, 
there  is  no  guarantee  that  the  hope  will  be  realized. 

The  iteration  Steps  4  through  7  are  now  repeated  until  convergence 
is  obtained. 

3.2  INTERPOLATION  AND  EXTRAPOLATION  SCHEMES 

During  the  course  of  the  iterative  procedure,  the  free  surface 
location  is  adjusted  by  satisfying  the  tangency  conditions  on  the  free 
surface  (Equation  2.6).  It  is  necessary,  subsequently,  to  regenerate  the 
segment  of  the  finite  element  mesh  adjacent  to  the  free  surface.  The  mid¬ 
side  nodes  on  the  free  surface  are  relocated  in  such  a  way  that  the  free 
surface  remains  smooth.  Those  nodes  lying  along  a  streamline  are  naturally 
and  smoothly  relocated  during  relocation  of  a  streamline  by  integration 
from  the  separation  point.  However,  an  interpolation  scheme  is  needed  to 
relocate  nodes  between  streamlines.  Because  of  the  symmetry  of  the  physical 
problem  for  the  pure  drag  case,  a  two-dimensional  interpolation  scheme  is 
sufficient.  However,  for  the  lifting  flow  case,  a  three-dimensional  scheme 
is  necessary. 


3.2.1  For  the  Pure  Drag  Case 

Referring  to  Figure  2.4,  one  sees  that,  if  the  movements  of  the  node 
points  on  the  free  surface  are  restricted  such  that  the  corner  nodes  always 
lie  on  planes  of  constant  x  ,  a  two-dimensional  interpolation  can  be  used 
to  locate  the  y  and  z  coordinates  of  mid-side  nodes  at  any  cross  section 
of  constant  x  given  that  the  location  of  the  corner  nodes  is  known. 

Because  of  the  symmetry  of  the  physical  problem  with  respect  to  the  y  and 
z  axes  as  shown  in  Figure  3.1,  a  periodic  cubic  spline  interpolation 
expressed  in  terms  of  the  circular  coordinates  0  and  r  is  chosen  with 
the  boundary  condition  r I *  r^0=2iT  *  However>  an  algorithm  for 

unequally  spaced  knots  (nodes)  is  required.  Since  the  theory  of  spline 
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Interpolation  is  well  documented,  readers  are  referred,  e.g.,  to  de  Boor 
(1978)  for  a  detailed  treatment  on  the  subject.  However,  the  algorithm 
used  in  the  present  study  is  taken  from  Spath  (1973,  Ch.  4). 

3.2.2  For  the  Lifting  Flow  Case 

In  order  that  the  effect  of  mid-side  node  locations  at  separation  as 
described  in  Section  2. 2. 3. 3  be  the  same  in  the  lifting  flow  at  all  itera¬ 
tions  as  in  the  pure  drag  case,  it  is  necessary  to  have  the  mid-side  node 
# 14  (Figure  2.6)  at  the  1/4  distance  from  the  separation  point  S  in  a 
transformed  x'  direction  which  is  perpendicular  to  the  plate  as  shown  in 
Figure  3.2.  Consequently,  it  is  desirable  that  the  movements  of  all  of 
the  node  points  on  the  free  surface  be  restricted  such  that  the  corner  nodes 
always  lie  on  planes  of  constant  x'  as  illustrated  in  Figure  3.2.  However, 
then  the  nodes  at  the  "end"  of  the  new  (relocated)  free  surface  may  fall 
short  of  or  over  the  centerline  ((^  )  so  that  extrapolation  or  interpolation, 
respectively,  is  necessary  to  align  the  corner  nodes  with  the  centerline. 
Unlike  the  pure  drag  case,  a  general  three  dimensional  algorithm  for  unevenly 
spaced  data  is  needed  here. 

Interpolation  with  irregularly  spaced  data  has  long  been  a  subject  of 
interest  for  numerical  analysts.  Among  the  methods  used,  one  can  find 
methods  based  on  local  bivariate  interpolations  (e.g.  Akima,  1978)  as  well 
as  methods  with  multivariate  smoothing  splines  (e.g.  Fulker,  1975).  These 
methods  in  general,  are  very  expensive  and  relatively  difficult  to  use. 
Following  the  concepts  inherent  in  the  isoparametric  finite  elements  as 
used  in  the  present  study,  a  "local",  general,  three-dimensional  Interpola¬ 
tion  and  extrapolation  scheme  by  means  of  isoparametric  mapping  has  been 
developed.  This  new  scheme  is  found  to  be  applicable  in  two,  as  well  as  in 
three,  dimensions. 

The  two-dimensional  result  is  obtained  first.  Consider  the  isopara¬ 
metric  transformation  given  in  Equation  2.26  and  the  shape  functions  given 
in  Equation  2.27  (cf.,  Equation  2.37  for  the  two  dimensional  shape  functions) 
and  Imagine  that  the  z  axis  in  the  physical  coordinates  is  identical  with 
the  5  axis  in  the  local  coordinates.  Then,  the  transformation  degenerates 
to  a  two-dimensional  problem.  Referring  to  Figure  3.3  and  mapping  the  side 
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=  1,  one  has  the  following 


x  *  x5  +  (xg  -  Xg)  -^  +  (x3  +  Xg  -  2Xg)  , 


y  -  y5  +  <y8  -  y3>  t  +  (y 3  +  ys  '  2y5}  2  » 


(3.1) 


-1  <  n  <  1 


where  x  and  y  are  the  coordinates  in  physical  space  for  any  point  along 
the  side  Imagine  the  side  dsoParametric 

finite  element  is  arranged  so  as  to  coincide  with  the  data  points  (i.e., 
the  known  node  points)  as  indicated  in  Figure  3.3(c).  Then,  one  can  use 
Equation  3.1  to  interpolate  the  value  of  y  (at  the  point  indicated  by 
the  cross,  for  example),  say,  for  a  known  value  of  x  between  x  and  x  , 
first,  by  solving  for  the  value  of  n  as 


where 


2 

axn  +  a2r'  +  a3  m  0 

al  =  ^x3  +  xg  “  2x5)/2  , 
a2  *  ^x8  “  x3^/2  »  and 


(3.2) 


and  second,  by  evaluating  y  from  the  second  equation  of  Equations  3.1. 

For  interpolation  between  in  Figure  3.3(c),  one  can  do  the  inter¬ 

polation  from  the  left  (using  C^“(^“CD  )  as  well  as  from  the  right 
(using  - /j-^  )  and  a  weighted  mean  value  from  the  two  interpolated 

values  can  be  taken  as  the  final  solution.  From  numerical  experiments,  it 

4 

is  observed  that  this  weighted  interpolation  scheme  is  approximately  0(S  ) , 
i.e.,  fourth  order  accurate  with  respect  to  the  side  length.  However,  this 
scheme  suffers  the  restriction  that  node  //5  must  lie  in  the  middle  half  of 
the  side  (T)-(T)-(T)  (°r  for  uniqueness  in  the  solution  to  be 

guaranteed. 
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One  can  also  extrapolate  by  means  of  the  Isoparametric  mapping. 
Specifically,  given  the  coordinates  of  and  the  cross  point  in 

Figure  3.3(c),  one  can  obtain  yg  provided  that  Xg  is  known  at  the 
extreme  point  #8.  As  soon  as  the  value  of  n  is  obtained  from  Equation  3.2, 
yg  can  be  found  from  the  following: 


y8 


y  +  (n2  -  l)  y5  +  ^  (i  -  n)y3 
§  (l  +  n) 


(3.3) 


where  y  and  q  are  values  for  the  intermdiate  point  x. 

Should  the  data  points  not  lie  on  the  same  z-plane,  one  has  to  add  the 
third  equation  3.1  for  z  ,  namely. 


z5  +  (Zg  -  z3>  f  +  (z3  +  z8  -  2z5)  5  , 


-l  <  n  <  l 


(3.4) 


where  z  is  the  z-coordinate  for  the  intermediate  point  at  n  .  By  means 
of  this  scheme,  one  can  do  interpolation  or  extrapolation  in  three  dimen¬ 
sions  with  little  computational  effort  and  any  one  of  the  x  ,  y  or  z 
coordinates  can  be  taken  as  the  independent  variable.  However,  this  is 
only  a  local  interpolation  scheme  such  that  any  changes  in  data  beyond  the 
two  points  on  either  side  of  the  interpolated  value  do  not  affect  the  final 
solution. 


3.3  WEIGHTING  SCHEME  AT  THE  FREE  SURFACE 

As  mentioned  in  Section  3.1,  the  separation  edge  represents  a  boundary 
singularity  (Crank  and  Furzeland,  1978)  at  which  two  types  of  boundary  con¬ 
ditions  meet,  namely,  the  no-flux  condition  on  the  solid  plate  and  the  pre¬ 
scribed  velocity  potential  <j)  from  the  free  surface.  Crank  and  Furzeland 
(1978)  discuss  the  nature  of  this  mixed  boundary  condition  singularity.  Such 
singularities  are  found  in  a  wide  variety  of  problems  and  various  authors 
have  used  a  variety  of  ways  to  combat  the  difficulties  introduced  by  the 
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singularity.  Of  particular  value  are  analytical  forms  which  allow  the 
singular  nature  of  the  problem  to  be  analytically  represented  within  the 
context  of  the  numerical  scheme.  For  two-dimensional  problems  with  regu¬ 
lar  boundaries  (boundaries  which  are  not  curved) ,  significant  success  has 
been  achieved  with  such  analytical  patches.  On  the  other  hand,  the  present 
cavity  flow  problem  has  both  curved  boundaries  and  irregularly-sized  finite 
elements  in  three  dimensions.  It  is  not  clear  that  an  analytical  patch  of 
the  nature  required  could  be  developed  for  three-dimensional  physical  space 
or  perhaps  in  the  isoparametric  space.  In  any  case,  an  alternative  approach 
used  by  many  authors  is  grid  refinement  near  the  singularity.  In  the  end, 
because  the  error  in  the  solution  for  velocity  potential  <J>  at  the  singu¬ 
larity  present  here  was  shown  by  Crank  and  Furzeland  (1978)  to  be  a  function 
of  nodal  distance,  it  was  decided  to  use  mesh  refinement  for  the  present 
solution.  Although  Crank  and  Furzeland  (1978)  demonstrated  specifically 
the  behavior  in  the  neighborhood  of  a  singularity  for  a  finite  difference 
solution  similar  behavior  occurs  for  finite  elements  (e.g.,  see  Strang  and 
Fix,  1973,  Sec.  8.4).  It  was  expected,  then,  that  by  reduction  of  the  finite 
element  size  in  the  neighborhood  of  the  separation  point  the  error  which  is 
manifested  in  the  present  problem  by  failure  to  satisfy  the  tangent-separa¬ 
tion  condition  exactly  could  be  systematically  reduced. 

The  separation  condition  requires  that  the  velocity  normal  to  the 
surface  of  a  body  at  separation  be  zero  and  that  the  flow  be  wholly  tangen¬ 
tial.  With  the  present  error  due  to  the  mixed  boundary  condition  singularity, 
the  normal  velocity  at  different  points  along  the  separation  edge  is  non-zero 
and  variable.  As  shown  in  Section  2.4,  the  free  surface  is  essentially  a 
characteristic  surface.  From  Section  2. 2. 3. 3  it  is  known  that  the  stream¬ 
lines  leave  the  separation  surface  tangential  to  the  solid  surface.  If  then 
the  separation  condition  is  not  satisfied  exactly,  the  governing  first-order- 
partial-differential-equation  for  the  free  surface  is  not  uniformly  satisfied 
at  the  separation  points.  This  occurs  because,  while  the  characteristic  lines 
(or  streamlines)  are  tangent  to  the  body,  they  are  not  tangent  to  the  velocity 
vector  there,  since  the  normal  velocity  at  the  separation  point  is  non-zero. 

If  the  separation  condition  is  not  uniformly  satisfied  at  each  point  along 
the  separation  edge,  the  streamlines  from  the  separation  edge  each  represent 
different  physical  flows.  Accordingly,  one  might  expect  progressive  divergence 
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of  the  iterative  solution.  Interestingly,  Larock  and  Taylor  (1976),  who 
made  no  effort  to  impose  a  tangent-separation  condition  for  their  jet  flow 
from  an  orifice,  found  that  their  solution  tended  to  drift  and  that,  under 
circumstances  in  which  the  flow  was  more  highly  three  dimensional,  conver¬ 
gence  was  either  slowed  or  not  obtained  at  all. 

It  became  necessary,  then,  to  devise  a  scheme  whereby  the  velocity 

potential  0  at  the  separation  edge  could  be  adjusted  to  control  satis- 
s 

faction  of  the  tangent-separation  condition.  This  led  to  the  development 
of  the  weighting  scheme  described  below.  From  the  computed  velocities  normal 
to  the  plate  at  various  separation  points,  one  can  derive  a  set  of  weights, 
one  for  each  separation  point  and  streamline,  which  adjust  the  velocity 
potential  at  separation  such  that  the  normal  velocity  at  each  separation 
point  is  zero  within  a  certain  limit.  The  limit  was  usually  set  as  follows: 


Max 


Velocity  normal  to  the 
plate  at  separation 


<  10 


-2 


(3.5) 


After  the  value  of  the  separation  velocity  potential  0g  is  determined, 

values  for  the  new  velocity  potential  along  the  appropriate  streamline  are 

apportioned  to  the  node  points  of  the  streamline  such  that  A0/As  constant. 

Since  <|)  equals  0  at  separation  and  0  is  zero  at  the  center  line  of 
s 

the  cavity,  the  free  surface  velocity  on  a  streamline  becomes  q  which  is 

s 

equal  to  <pg  divided  by  the  length  of  the  streamline  in  question.  Then 

the  new  set  of  free-surface  velocity  potentials  0^  are  used  as  the 

Dirichlet  boundary  conditions  on  the  free  surface  for  the  next  iteration; 

however,  since  0  is  determined  to  make  the  normal  velocities  zero,  q 
s  s 

is  not  necessarily  the  same  on  each  streamline.  This  is  the  genesis  of  the 
main  difficulties  encountered  herein. 

This  weighting  scheme  produced  an  absolutely  convergent  iterative 
process.  Other  schemes  which  we  employed  gave  unstable  or  divergent 
iterative  processes  except  in  those  cases  in  which  the  flow  was  two-dimensional, 
i.e.,  axisymmetric  or  plane  flows.  The  drifting  of  the  streamlines  reported 
by  Larock  and  Taylor  (1976)  did  not  occur.  Detailed  results  are  discussed  in 
Sections  4.2,  5  and  6;  however,  the  general  procedure  in  solving  problems 

using  the  weighting  scheme  is  as  follows.  To  start,  the  upstream  flow  is 
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prescribed  and  held  constant  throughout  the  procedure.  With  the  assumed 

location  of  the  free  surface,  the  free  surface  velocity  can  be  estimated 

from  a  continuity  argument  applied  to  the  upstream  face  of  the  computation 

region  and  the  center  line  of  the  cavity.  Weights  are  equal  to  unity  in  the 

first  iteration;  therefore,  the  velocity  potential  <f>s  on  the  free  surface 

at  separation  equals  the  length  of  a  streamline  times  the  estimated  constant 

velocity  on  the  free  surface.  The  problem  for  $  is  then  solved,  and  the 

velocities  at  separation  are  computed.  Based  on  the  magnitude  of  the  velocity 

normal  to  the  plate,  the  weights  are  adjusted  such  that,  when  they  are 

multiplied  with  the  original  velocity  potential  <f>  at  separation,  the  sepa- 

'  *  s 

ration  condition  should  be  satisfied.  For  example,  if  there  is  a  positive 
normal  velocity  in  a  pure  drag  flow,  such  as  that  depicted  in  Fig.  2.4,  the 
velocity  potential  $  will  be  increased  to  reduce  this  positive  normal 
velocity.  These  new  values  of  the  velocity  potential  (j)^  at  separation  are 
then  prorated  to  the  nodes  on  each  streamline  proportional  to  the  distance 
along  the  streamline  as  indicated  above.  The  problem  is  now  solved  again. 

This  procedure  is  repeated  until  the  separation  criterion  is  satisfied,  and 
the  free  surface  is  then  moved.  With  the  new  location  of  the  free  surface, 
the  resulting  new  problem  is  solved  again  (as  above)  until  convergence  is 
obtained.  Because  continuity  is  satisfied  after  one  or  two  iterations,  the 
solution  is  considered  to  have  converged  when  the  stream  surface  stops 
moving.  This  is  normally  taken  to  occur  when  no  change  occur  in  the  fourth 
significant  figure  for  the  physical  location  of  any  free  service  point  on 
the  cavity  center  line,  i.e.,  at  the  downstream  end  of  the  computational 
area,  over  several  iterations.  During  the  iteration  procedure,  it  was  often 
found  to  be  more  economical  to  adjust  the  weights  manually  between  each 
iteration.  This  enables  the  user  to  lead  the  iterative  process  to  more 
rapid  convergence  and  reduce  the  cost  of  the  computation.  In  the  cases 
studied  here,  the  circular  arc  and  elliptic  pure  drag  cases  were  iterated 
with  user  itervention  until  the  first  move  of  the  free  surface.  Automatic 
correction  was  used  thereafter.  Lifting  flow  cases  were  iterated  manually. 

See  Ko  and  Street  (1979) . 


4.  IMPLEMENTATIONS  AND  TESTS 

4.1  SPECIFIC  IMPLEMENTATIONS 

In  Section  2.2.3,  the  implementation  of  the  boundary  conditions  was 
addressed.  It  was  found  that  the  condition  that  the  free  surface  be  tangent 
to  the  surface  of  the  foil  at  the  separation  point  is  satisfied  by  the 
quarter  point  location  of  the  first  mid-node.  This  is  true  asymptotically 
for  the  particular  set  of  shape  functions  chosen.  However,  another  method 
was  considered,  namely,  patching  in  an  analytic  solution  for  the  free  sur¬ 
face  shape  (at  the  separation  point)  similar  to  the  one  used  by  Mogel  and 
Street  (1974).  In  two-dimensional  flow,  the  form  of  the  free-surface  near 
the  separation  point  varies  with  the  two-thirds  power  of  the  distance  away 
from  the  point,  i.e.. 


y 


2/3 


(4.1) 


with  equality  given  by 


y  -  Cx2/3  (4.2) 

where  C  ,  the  constant  of  proportionality,  is  determined  solely  by  the 
velocity-components  at  the  free  surface  at  a  horizontal  distance  x  from 
the  separation  point.  Assuming  that  the  flow  near  the  separation  point  is 
essentially  two-dimensional  (locally),  the  relationship  given  in  Equation 
4.1  is  used  in  the  three-dimensional  case  with  y  being  measured  perpen¬ 
dicular  to  the  edge  of  the  foil  and  in  a  direction  tangent  to  the  foil. 

This  analytic  patching  algorithm  was  implemented  and  tested  with  the  pure 
drag  flow  cases.  It  was  found,  oftentimes,  that  analytic  patching  caused 
initially  a  very  sharp  curvature  of  the  free  surface  such  that  a  smooth 
continuation  of  the  rest  of  the  free  surface  was  impossible.  When  converg¬ 
ence  was  obtained,  the  patch  had  no  significant  effect  on  the  results  com¬ 
pared  to  cases  without  the  patch.  Therefore,  the  patch  technique  was 
dropped  from  the  program  as  an  option  for  free  surface  moving  routine. 
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During  Che  iterative  procedure,  the  free  surface  Is  relocated  based 
on  the  tangency  of  the  velocity  vector  to  the  free  surface.  This  is 
achieved  by  satisfying  Equation  2.6  in  the  form 


(4.3) 


Equations  4.3  can  be  solved  by  direct  integration,  with  the  intermediate 
values  of  u  ,  v  and  w  being  interpolated  linearly  between  node  points, 
which  is  consistent  with  the  quadratic  finite  element  approximation.  A 
fourth-order  Runge-Kutta  procedure  is  used.  However,  at  the  first  element 
next  to  the  separation  edge,  the  inverse  of  the  Equation  4.3  has  to  be  used 
to  avoid  infinite  slope  at  the  separation,  namely, 


(4.4) 


By  means  of  Equations  4.4,  one  can  also  specify  the  zero  slope  conditions, 
namely , 

dx  dx  n  ,,  -v 

0  (4.5) 

dy  dz 

at  the  separation  point.  However,  this  did  not  seem  to  be  crucial  for  the 
iterative  procedure  to  converge. 

4.2  VERIFICATION  OF  ALGORITHMS 

In  order  to  test  the  performance  of  the  selected  algorithms  the 
present  procedure  has  been  tested  against  the  best  data  available.  Although 
there  is  not  exact  solution  for  the  present  three-dimensional  cases  with 
which  to  compare,  it  is  hoped  that  the  individual  tests  performed  give 
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insight  into  the  overall  performance  of  the  procedure.  Two  types  of  tests 
were  undertaken.  First,  two-dimensional  flows  were  tested,  viz.,  a  plane 
cavity  flow  and  an  axisymmetric  flow  were  solved  by  use  of  the  three- 
dimensional  program.  Second,  the  free  surface  shift  algorithm  was  tested 
alone  for  a  fully  three-dimensional  flow. 

4.2.1  Two-Dimensional  Geometries 

Because  the  flow-field  has  some  basic  symmetry,  no  variations  are  to 
be  expected  with  depth  (that  is  across  the  flow)  in  the  plane  flow  or  in 
azimuth  for  the  axisymmetric  flow.  Accordingly,  the  weighting  scheme  of 
Section  3.3  was  not  employed.  The  normal  velocity  at  separation  approaches 
zero  then  only  as  a  consequence  of  the  natural  iteration  process  and  is  not 
controlled  or  otherwise  constrained  except  by  the  adjustment  of  the  free 
surface  velocity  to  yield  balanced  inflow  and  outflow.  The  velocity  poten¬ 
tial  on  the  free  surface  is  prescribed  then  according  to  Equation  2.5  for 
every  iteration. 

One  two-dimensional  geometry  was  selected  from  among  those  reported 
in  Mogel  and  Street  (1974)  as  an  example  of  the  behavior  of  the  three- 
dimensional  program  for  a  two-dimensional  flow.  As  shown  in  Figure  4.1  a 
flow  was  constructed  which  has  a  typical  cavity  configuration  in  the  plane 
of  flow  and  is  one  element  deep.  If  the  program  is  operating  correctly, 
there  are  no  driving  forces  to  cause  the  front  and  back  faces  of  the  element 
with  respect  to  depth  to  have  different  values  of  the  unknown  potential  or 
velocities.  Accordingly,  the  correct  program  will  lead  to  a  convergent 
solution  in  which  the  values  at  the  corresponding  points  of  the  front  and 
back  faces  of  any  element  are  essentially  equal  and  the  flow  remains  two- 
dimensional.  When  tested  the  present  procedure  yielded  a  converged  solution 
in  23  iterations.  The  criterion  for  convergence  was  defined  as 

£  <  10"2  (4.6) 

with  e  being  the  maximum  movement  of  any  node-point  on  the  free  surface 
between  two  consecutive  iterations.  Figure  4.2  shows  the  sequence  of  the 
movement  of  the  free  surface  in  the  present  case.  Interestingly  the  three- 
dimensional,  finite  element  code  is  found  to  be  no  more  expensive  for  this 
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case  than  the  two-dimensional  method  used  by  Mogel  and  Street  (1974)  for 

a  solution  of  similar  accuracy.  However,  in  their  paper  it  was  possible 

to  prescribe  exactly  that  the  x-component  of  the  fluid  velocity  vector  at 

separation  from  the  plate  is  identically  equal  to  zero.  In  the  present 

case  the  shape  function  used  in  the  finite  element  method  leads  to  an 

_2 

x-component  of  the  velocity  at  separation  of  the  order  of  10  times  the 
maximum  velocity  on  the  free  surface. 

Working  in  stream  function-potential  space  Brennen  (1969)  solved  the 
axisymmetric  cavity  flow  past  a  circular  disk  in  a  circular  water  tunnel. 

He  used  a  Riabouchinsky  image  model  and  a  finite  difference  numerical 
technique.  His  numerical  results  were  verified  by  a  series  of  experiments. 
It  is,  therefore,  a  good  test  of  the  present  procedure  to  compare  the 
finite  element  result  with  Brennen* s  results  for  a  typical  case.  Actually, 
three  cases  with  different  geometric  dimensions  were  chosen  for  the  tests. 
Converged  solutions  for  all  three  cases  were  achieved  with  no  weighting 
(cf . ,  Section  3.3)  . 

The  results  of  these  tests  are  shown  in  Figure  4.3  together  with  the 
curves  presented  by  Brennen  (1969).  It  is  clear  that  the  finite  element 
and  Brennen  techniques  are  in  precise  quantitative  agreement.  In  addition 
this  test  shows  that  if,  as  a  consequence  of  flow  geometry  the  discretiza¬ 
tion  through  the  finite  element  prescription  occurs  unfiormly  across  the 
flow,  then  errors  in  flow  separation  velocity  likewise  occur  uniformly  and 
the  weighting  scheme  described  in  Section  3.3  is  not  necessary  in  order  for 
the  iteration  procedure  to  be  stable;  however,  it  leads  to  the  same  results 
as  obtained  here. 

4.2.2  Three-Dimensional  Geometry 

The  analytic  solution  for  a  uniform  flow  past  a  three-dimensional 
Rankine  Oval  is  well-known.  In  particular,  for  a  particular  prescription 
of  the  oval  geometry,  i.e.,  length  and  width,  one  knows  the  exact  surface 
geometry  of  the  oval  and  the  entire  velocity  field  surrounding  it.  It  is 
possible  then  to  use  this  information  to  test  the  present  algorithm  for 
moving  the  free  surface  without  being  involved  with  the  discretization, 
truncation  or  computation  errors  inherent  in  solving  for  the  potential  and 


velocity  fields  in  a  given  case.  We  begin  by  generating  a  regular,  three- 
dimensional,  geometric  mesh  for  a  typical  cavity  flow  past  a  disk  in  a 
water  tunnel,  for  example,  as  shown  by  the  solid  lines  in  Figure  4.4. 

Then  the  three-dimensional  Rankine  Oval  is  placed  in  the  cavity  such  that 
the  edge  of  the  disk  coincides  with  the  analytic  boundary  of  the  oval  at 
some  convenient  point.  At  each  iteration  the  velocity  field  is  thus  com¬ 
puted  from  the  analytic  solution  for  the  Rankine  Oval  and  the  free  surface 
is  moved  from  its  original  position  in  accordance  with  the  analytic  velocities 
for  the  Oval  and  the  free  surface  moving  algorithm  as  shown  in  Figure  4.4. 

The  moved  free  surface  converges  to  the  analytic  shape  of  the  Rankine  Oval. 

The  accuracy  obtained  with  a  reasonable  number  of  iterations  is  of  the  order 
of  four  to  six  significant  figures. 

From  the  results  of  the  above  tests,  it  is  clear  that  the  finite 
element  method  is  more  efficient  than  the  finite  difference  method  used  by 
Mogel  and  Street  (1974) .  From  the  tests  of  the  plane  and  axisymmetric  two- 
dimensional  flows  it  is  clear  that  the  iterative  procedure  is  absolutely 
convergent  in  those  cases  in  which  the  representation  of  the  separation 
velocity  at  the  start  of  the  free  surface  is  not  influenced  by  true  three- 
dimensional  effects.  Further,  the  Rankine  Oval  test  shows  that  if  the 
velocity  field  is  accurately  represented  the  iterative  procedure  for  moving 
the  free  stream  surface  is  rapidly  and  stably  convergent  to  any  desired 
level  of  accuracy.  We  will  see,  however,  in  Sections  5.1  and  5.2  that 
discretization  errors  and  the  variations  introduced  by  truly  three-dimensional 
flows  (those  about  strongly  elliptic  plates  and  lifting  flows)  can  cause 
either  an  unstable  iteration  or  convergence  to  a  solution  in  which  the  free 
stream  surface  pressure  is  not  constant  as  one  moves  across  or  transverse 
to  the  streamlines  of  the  flow. 
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5.  CAVITY  FLOW  RESULTS 

One  of  the  objectives  in  developing  the  present  model  is  to  gain  the 
ability  to  simulate  nonlinear,  fully  cavitating  flows  about  hydrofoils.  As 
mentioned  in  Section  1,  there  is  not  a  fully  nonlinear,  three-dimensional 
model  that  can  be  used  for  design  purposes.  Consequently,  design  is  often 
based  on  water  tunnel  experimental  data.  However,  the  data  collected  from 
water  tunnel  experiments  is  strongly  related  to  the  physical  configuration 
of  the  experimental  setup.  The  present  model  makes  it  possible  to  examine 
the  effects  of  different  physical  configurations.  Thus,  one  has  direct 
insight  into  the  interpretation  of  experimental  results. 

In  this  section  results  are  presented  for  three  flow  geometries.  In 
Section  5.1  pure  drag  flows  past  circular  and  elliptic  plates  in  rectangular 
water  tunnels  are  described.  A  lifting  flow  consisting  of  a  circular  plate 
at  an  angle  of  attack  of  60°  in  a  rectangular  water  tunnel  is  given  in 
Section  5.2 

Before  examination  of  these  flows  it  is  crucial  to  understand  the 
accuracy  of  presented  solutions,  which  will  be  discussed  in  detail  in 
Section  6.  Because  of  the  separation  edge  singularity  reviewed  in  Section 
3.3,  the  ability  of  the  numerical  solution  to  maintain  a  constant  pressure 
(or  velocity)  over  the  entire  cavity  decreases  as  the  three  dimensionality 
of  the  free  surface  increases .  The  pressure  (velocity)  is  closely  constant 
streamwise  or  along  each  streamline;  however,  as  one  moves  across  the  flow 
transverse  to  the  streamlines  the  pressure  (velocity)  changes.  The  causes 
and  potential  cures  for  this  effective  error  in  the  solution  are  covered  in 
Section  6. 

The  effective  error  is,  thus,  negligible  for  all  of  the  circular 
plate  pure  drag  cases  in  which  the  cavity  cross  section  is  close  to  being 
circular.  These  solutions  are  essentially  exact.  The  total  difference 
between  free  stream  velocities  in  the  elliptic  plate  cases  is  of  the  order 
of  five  percent  (10  percent  for  the  pressure  coefficient).  However,  for 
the  lifting  flow,  the  total  free  stream  velocity  difference  between  the 
leading  and  trailing  edge  streamlines  is  about  10  percent  (21  percent  for 
the  pressure)  which  is  unacceptable. 
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5.1  RESULTS  FOR  THREE-DIMENSIONAL,  PURE  DRAG  CASES 


5.1.1  Circular  Plates  in  Square  Water  Tunnels 

A  circular  plate  with  a  two  unit  radius  is  used  with  three  different 
square  water  tunnels  with  widths  of  10,  15  and  20  units.  This  configuration 
was  chosen  because  it  is  the  simplest  variation  from  an  axisymmetric  flow 
with  a  circular  plate  in  a  circular  water  tunnel  (as  solved  by  Brennen,  1969; 
see  Figure  4.3).  The  circular  plate  is  placed  at  a  distance  of  10  units 
from  the  upstream  boundary  at  which  a  uniform  flow  is  assumed  (i.e. ,  F  *  10; 
cf.,  Figure  2.4).  Due  to  the  symmetry  of  the  setup  as  mentioned  in  Section 

2.1,  only  a  quarter  of  the  flow  field  needs  to  be  solved.  For  the  water 
tunnels  with  widths  of  10,  15  and  20  units,  the  half  width  of  the  square 
tunnel  W  takes  on  values  of  5,  7.5  and  10  respectively.  Two  different 
half-lengths  L  for  the  cavity  are  used,  namely,  5  and  10  units.  The 
effect  of  changing  the  cavity  length  will  be  discussed  later  in  this 
section.  For  the  case  with  a  water  tunnel  width  of  10  units,  the  effect  of 
grid  refinement  was  also  examined  (See  Section  6.2).  The  apparent  effect 
on  the  results  was  small,  but  some  conclusions  are  drawn  later  in  Section 

6.2. 

The  results  can  be  viewed  from  several  perspectives,  namely,  the 
cavity's  gross  shape,  the  wall  effects,  the  "stability  of  pressure  coeffici¬ 
ent"  observed  by  Wu,  et  al.  (1971)  and  Mogel  and  Street  (1974),  and  the 
accuracy  of  the  solution.  These  are  treated  in  turn  in  the  following 
paragraphs . 

Interestingly,  in  all  cases  examined  here,  the  cavity  remains  rela¬ 
tively  circular  in  spite  of  the  influence  of  the  square  water  tunnel  (cf.. 
Figures  5.1  and  5.2).  (For  reference  the  lines  are  circular  arcs  in  these 
figures.)  However,  the  wall  effect  is  substantial. 

Normally,  one  views  results  at  constant  cavitation  number  a  and 
with  various  tunnel  width  to  plate  radius  ratios.  Then,  a  decrease  in  tunnel 
size  produces  an  increase  in  cavity  length  for  a  given  o  as  the  flow  is 
driven  toward  choked  flow.  However,  it  is  more  convenient  with  a  physical- 
space  numerical-scheme  to  hold  cavity  length  2L  constant  and  let  cavitation 
number  0  vary.  With  cavity  length  held  constant,  decreasing  the  tunnel 
size  tends  to  increase  0  .  This  results  also  in  a  decrease  in  the  B/P  ratio 
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(cf.,  Figure  2.4).  The  corresponding  increase  in  the  velocity  on  the  free 
surface  is  observed  in  Table  5.1.  Figure  5.3  shows  two  sets  of  curves  for 
the  square  water  tunnel  with  two  different  cavity  lengths,  namely,  L  *  5 
and  10  units  with  L/P  *  2.5  and  5.0,  respectively.  When  W/P  is  small, 
a  large  a  is  needed  to  maintain  as  a  short  cavity;  then  the  actual  cavity 
length  becomes  relatively  unimportant  and  the  cavitation  numbers  for  the 
square  tunnel  should  be  close  together  (as  shown  in  Figure  5.3)  because  the 
L/P  ratios  studied  are  different  only  by  a  ratio  of  two.  This  is  also  true 
for  the  ratio  B/P. 

From  Figures  5.1  and  5.2,  one  might  conclude  that  the  comer  of  the 
square  water  tunnel  had  little  effect  on  the  cavity  shape.  However,  the 
blockage  effect  is  very  different  from  a  similar  axisymmetric  case.  Figure 
5.3  demonstrates  these  trends.  For  large  W/P  ,  the  wall  effects  are  insig¬ 
nificant  and  the  blockage  effect  of  a  circular  tunnel  is  no  different  than 
that  of  the  square  tunnel.  Therefore,  the  cavitation  number  a  from  the 
square  tunnel  coincides  with  the  similar  axisymmetric  case  for  both  values 
of  L/P  tested  here.  However,  when  the  width  of  the  square  tunnel  (the 
radius  in  the  circular  case)  decreases,  the  axisymmetric  flows  feel  the  wall 
effect  much  faster  than  the  similar  square  tunnel  flows.  For  the  axisymmetric 
flows,  the  cavitation  numbers  a  converge  for  the  two  L/P  ratios  when 
W/P  decreases,  as  observed  for  the  square  tunnel  flows. 

Wu  et  al.  (1971)  and  Mogel  and  Street  (1974)  observed  that  the  coeffici¬ 
ent  C*  •  C^d  +  cr)  ^  is  nearly  constant  for  many  flow  configurations  in 
which  W/P  is  varied.  Table  5.1  shows  that  C*  is  close  to  constant  in  this 
three  dimensional  flow  as  well.  However,  while  C*  appears  virtually  inde¬ 
pendent  of  a  ,  C*  does  increase  by  about  5%  when  W/P  is  reduced  by  50%. 

Table  5.1  summarizes  the  results  for  a  particular  set  of  circular- 

plate-in-square-tunnel  cases.  For  each  case,  convergence  was  certainly 

achieved  in  25  iterations;  actually  no  significant  changes  are  observed  after 

15-20  iterations.  Convergence  means  that  the  velocity  components  normal  to 

-2 

the  plate  at  separation  were  less  than  about  10  times  the  free  stream 
velocity  and  no  changes  were  observed  in  the  fourth  significant  figure  of  the 
cavity  streamline  end  points  over  several  iterations.  Contrary  to  the  experi¬ 
ence  of  Larock  and  Taylor  (1976),  no  accumulated  drift  in  the  streamlines  is 
observed  in  the  present  study. 
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5.1.2  Elliptic  Plates  in  Rectangular  Water  Tunnels 

The  elliptic  plate  used  by  Jiang  and  Leehey  (1977)  in  their  lifting 
flow  study  was  chosen  for  this  pure  drag  flow  study.  The  basic  configura¬ 
tion  consisted  of  an  elliptic  plate  with  2.35  unit  by  7.75  unit  semiaxes 
and  a  40  unit  by  20  unit  rectangular  water  tunnel.  With  reference  to 
Figure  2.4,  the  plate  is  placed  at  a  distance  of  20  units  from  the  upstream 
boundary  at  which  a  uniform  flow  is  assumed,  e.g.,  F  *  20  units.  In  order 
to  examine  the  wall  effect,  a  water  tunnel  with  80  unit  by  40  unit  was  also 
used.  Two  different  cavity  half-lengths,  namely,  10  units  and  20  units, 
were  used  to  examine  the  change  of  the  cavitation  number  with  respect  to 
the  tunnel  geometry.  The  parameters  of  the  studied  geometry  and  the  gross 
results  are  given  in  Table  5.2. 

In  spite  of  the  influence  of  the  rectangular  water  tunnel  shape,  the 
cavity  remains  essentially  elliptic  (Figure  5.4).  Again,  the  wall  effect 
is  substantial.  As  indicated  in  Figure  5.4,  as  the  ratio  of  the  tunnel  half 
width  to  the  short  semiaxis  of  the  plate  W/P  decreases  from  6.15  to  3.08, 
the  cavity  Is  pushed  in  toward  the  center  line,  and  there  is  an  increase  in 
the  velocity  on  the  free  surface,  i.e.,  a  increases  (cases  El  and  E3, 

Table  5.2).  However,  the  wall  effect  on  the  cavity  size  due  to  the  decrease 
In  the  ratio  of  D/Q  from  5.16  to  2.58  is  very  small. 

In  the  elliptic  cases  tested  here,  the  coefficient  C*  showed  a 
slight  dependence  on  tunnel  size,  but  not  on  o  ,  the  difference  being 
approximately  5 %  between  the  two  tunnel  configurations.  However,  this 
difference  may  be  due  to  the  error  in  the  constant  pressure  condition  on 
free  surface  (see  Section  6) . 

5.2  RESULTS  FOR  A  THREE-DIMENSIONAL  LIFTING  FLOW 

The  final  step  of  the  present  study  was  simulation  of  a  lifting  super- 
cavitating  flow  in  a  water  tunnel.  The  geometric  configuration  is  essentially 
the  same  as  the  pure  drag  flow  with  a  circular  plate  of  2  unit  radius  in  a 
square  water  tunnel  of  15  by  15  units  (Section  5.1.1).  The  center  of  the 
plate  is  placed  10  units  downstream  from  the  upstream  boundary  where  an 
uniform  velocity  U  is  prescribed.  The  plate  is  at  an  angle  of  attack  a 


of  60°.  The  cavity  half-length  is  5  units  from  the  center  of  the  plate. 
Figure  5.5  shows  the  FE  mesh  in  three  dimensions.  Figure  5.6  shows  the 
final  centerline  cavity  shape  after  5  shifts  of  the  free  surface.  The 
cross-section  of  the  cavity  is  found  to  be  only  a  slight  variation  from 
circular  (cf.,  the  dashed  line  in  Fig.  5.6).  The  stagnation  point  (where 
the  pressure  attained  its  maximum)  is  approximately  l/8th  of  the  diameter 
of  the  plate  from  the  leading  edge.  Figure  5.7  shows  the  pressure  distri¬ 
bution  along  the  centerline  of  the  plate.  Both  because  the  separation  edge 
represents  a  boundary  singularity  and  because  the  normal  method  of  computing 
the  flow  velocity  breaks  down  at  the  separation  edge  (see  Section  6) ,  the 
pressure  coefficient  on  the  plate  does  not  drop  off  to  zero  at  separation 
as  it  should.  The  effect  of  the  velocity  calculation  is  not  large  and  can 
be  corrected  for  if  the  boundary  singularity  effects  were  minimized. 

Figure  5.8  shows  a  contour  plot  of  the  pressure  distribution  on  the 
plate  for  the  lifting  flow.  Because  in  the  plot  straight  lines  are  used 
to  join  the  node  points  the  circular  edge  of  the  plate  appears  to  have 
comers.  The  intermediate  values  of  the  contours  are  interpolated  from 
those  at  the  nearby  node  points.  This  when  coupled  with  the  behavior  of 
the  element  interpolation  functions  and  the  contour  routine  interpolation 
functions  introduce  slight  wavy  variations  on  the  contours.  These  variations 
would  be  eliminated  if  more  node  points  were  used  on  the  plate.  The  cavita¬ 
tion  number  a  for  this  lifting  flow  is  0.53.  The  corresponding  drag  and 
lift  coefficients  C*  and  C*  are  0.63  and  0.37  respectively. 

The  lifting  flow  tested  here  converged  to  the  final  solution  after  5 
free  surface  relocations.  However,  it  took  approximately  30  iterations  to 
achieve  this  result.  It  was  clear  that  it  is  more  economical  to  manually 
adjust  the  "weights"  (cf.  Section  3.3)  between  iterations.  By  this,  human 
experience  can  be  employed  with  judgment  to  make  the  rate  of  convergence 
faster.  The  test  problem  has  303  elements  with  3137  node  points.  Each 
iteration,  including  setting  up  the  elemental  stiffness  matrices  for  all 
elements,  took  approximately  5-1/2  minutes  CPU  time  on  the  IBM  370/168 
computer.  The  core  requirement  was  1152  K  with  approximately  1500  tracks 
of  disk  space  on  IBM  3300  disks  for  temporary  storage. 
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6.  PERFORMANCE  OF  THE  SOLUTION  AND  ALTERNATIVE  FORMULATIONS 


The  computation  procedure  outlined  in  Sections  2  and  3  is  a  relatively 
straightforward  one.  The  program  structure  and  basic  finite  element  compu¬ 
tational  algorithms  were  adapted  from  an  existing  three-dimensional  ground- 
water  flow  code  for  free  surface  flows  developed  by  Dr.  C.  Taylor  at  the 
University  of  Wales,  Swansea.  The  crucial  choices  and  developments  focused 
on  the  iterative  method  used  herein  and  the  means  of  prescribing  and  dealing 
with  the  high  speed  free  surface,  which  has  a  different  boundary  condition 
than  that  found  in  groundwater  flow. 

Tests  to  verify  the  performance  of  the  computational  algorithms  and  to 
demonstrate  that  the  computational  method  reproduced,  exactly,  existing  axi- 
symmetric  flow  calculations  were  presented  in  Section  4.2.  The  capacity  of 
the  method  to  solve  fully  three-dimensional  cavity  flows  was  demonstrated 
in  Sections  5.1  and  5.2. 

It  is  now  appropriate  to  examine  the  nature  of  the  iterative  process 
that  was  used  in  this  work,  to  discuss  in  detail  the  accuracy  of  the  three- 
dimensional  solutions  presented  in  Section  5,  and  to  explore  alternative  formu¬ 
lations  which  might  remove  some  of  the  difficulties  encountered  with  the 
present  techniques.  These  three  areas  are  explored  in  Sections  6.1  through 
6.3. 


6.1  THE  NATURE  OF  THE  ITERATIVE  PROCESS 

In  order  to  separate  the  nonlinear  cavity  problem  into  linear  pieces 
it  is  necessary,  first,  to  establish  a  trial  location  of  the  cavity  free 
surface,  and  second,  to  prescribe  on  that  surface  the  velocity  potential  <J>  . 
When  the  solution  is  complete,  a  correct  and  exact  solution  must  have  a 
constant  velocity  along  each  streamline  of  the  free  surface,  and  the  velocity 
on  each  streamline  must  be  the  same  as  the  velocity  on  adjacent  streamlines. 
That  is,  there  will  be  a  single  constant  velocity  over  the  entire  free  sur¬ 
face.  In  addition,  the  actual  free  surface  must  be  tangent  to  the  solid 
plate  at  the  separation  edge.  While  this  can  be  mechanically  guaranteed, 
i.e.,  the  finite  element  method  forces  the  physical  surface  to  be  tangent  at 
the  body  at  separation  (see  Section  2. 2. 3. 3)  it  is  difficult  with  the  present 
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method  to  obtain  the  condition  on  the  velocity  that  it  be  tangent  to  the 
solid  body  at  separation.  Both  conditions  must  be  met  on  an  actual  free 
surface.  Accordingly,  in  formulating  the  problem  one  is  faced  with  two 
requirements:  first,  create  a  solution  in  which  the  velocity  is  constant 
over  the  entire  cavity;  second,  create  a  solution  in  which  the  velocity  is 
tangent,  not  only  to  the  free  surface  of  the  cavity,  but  also  to  the  solid 
body  at  separation. 

While  it  is  obvious  that  the  velocity  along  a  given  streamline  should 
be  constant  and,  hence,  that  the  potential  along  the  streamline  should  be 
specified  so  that  the  derivative  of  the  potential  with  regard  to  distance 
along  the  streamline  is  a  constant,  it  is  not  clear  that  one  can  simultane¬ 
ously  require  the  velocity  on  each  streamline  to  be  equal  to  the  same 
constant  value.  Thus,  in  this  work  two  alternative  procedures  for  specifying 
the  velocity  potential  were  examined.  Only  one  has  been  described  previously 
because  the  alternative  procedure,  namely,  that  used  by  Mogel  and  Street 
(1974)  and  (in  a  three-dimensional  flow  similar  to  the  present  one)  by 
Larock  and  Taylor  (1976) ,  led  to  unstable  solutions  for  the  present  cavity 
flow  problems.  This  unstable  method  consisted  of  adjusting  the  potential 
on  the  free  surface  such  that  the  velocity  was  constant  everywhere  on  the 
free  surface.  The  velocity  was  adjusted  so  as  to  promote  conservation  of 
mass  between  the  upstream  and  downstream  cross-sections  of  the  cavity  flow. 

No  attempt  was  made  to  control  the  tangent  separation  condition  or  to  force 
the  velocity  normal  to  the  plate  at  separation  to  be  zero.  The  alternative 
iteration  procedure,  which  was  the  one  adopted  here,  is  described  in  detail 
in  Section  3.3.  There,  while  the  velocity  was  kept  constant  on  each  indi¬ 
vidual  streamline,  it  was  allowed  to  vary  from  streamline  to  streamline 
because  the  requirement  was  Imposed  on  the  calculation  that  the  tangent 
separation  condition  be  satisfied.  This  led  to  an  exceedingly  stable  itera¬ 
tion.  However,  in  cases  of  highly  three-dimensional  flow,  that  is,  when 
the  cavity  shape  was  three-dimensional,  the  converged  solution  had  different 
velocities  on  adjacent  streamlines.  The  size  of  this  error  is  discussed  in 
the  next  section. 

The  question  arises:  Why  is  the  ostensibly  more  logical  iterative 
procedure  unstable?  An  examination  of  the  detailed  print-out  for  a  typical 
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convergent  simulation  by  Larock  and  Taylor  (1976)  [kindly  provided  by 
Professor  Larock]  shows  that  their  solution  does  not  satisfy  the  tangent 
separation  condition;  the  velocity  normal  to  the  orifice  of  their  jet  flow 
is  of  the  order  of  0.2  to  0.5  times  the  actual  velocity  at  separation,  and 
the  velocity  normal  to  the  plate  is  different  at  the  top  and  bottom  of  their 
orifice.  For  slightly  three-dimensional  flows  such  as  theirs  (and  indeed 
our  circular-plate  pure-drag  cases)  this  appears  to  be  of  no  consequence; 
however,  we  conjecture  that  failure  to  satisfy  the  tangent  separation  condi¬ 
tion  uniformly  may  be  a  root  cause  of  the  failure  of  the  lower  Froude-number 
cases  of  Larock  and  Taylor  (1976)  to  converge.  These  flows  are  more  highly 
three-dimensional  and  would  appear  to  suffer  from  the  same  errors  encountered 
in  the  present  work.  In  the  iteration  used  in  the  present  work  the  tangent 
separation  condition  is  satisfied  essentially  exactly  in  the  iteration  pro¬ 
cess  and  the  velocity  normal  to  a  solid  boundary  is  always  at  least  two 
orders  of  magnitude  smaller  than  the  average  velocity  at  separation. 

Larock  and  Taylor  (1976)  noted  that,  if  they  continued  to  iterate  their 
solution  and  shift  their  free  surface  as  indicated  by  their  iterative 
procedure,  a  number  of  small  shifts  tended  to  accumulate  to  a  non-negligible 
effect.  Thus,  their  solution  tended  to  drift.  If  one  uses  the  iterative 
procedure  developed  here  and  applies  it  to  a  flow  from  an  orifice  without  the 
effect  of  gravity  (a  problem  somewhat  similar  to  that  considered  by  Larock 
and  Taylor  at  high  Froude  numbers),  no  drift  of  the  solution  occurs.  Indeed, 
for  the  pure  drag  slightly  three-dimensional  cases  examined  here  the  itera¬ 
tive  procedure  described  in  Section  3.3  works  exceedingly  well  and  gives 
solutions  that  are  quantitatively  accurate.  The  present  method,  then,  may 
be  considered  perfectly  adequate  for  slightly  three-dimensional  flows. 
However,  the  problem  remains,  namely,  what  is  the  appropriate  formulation 
for  more  highly  three-dimensional  flows  in  which  the  competing  free  surface 
and  separation  boundary  conditions  really  must  be  satisfied  simultaneous? 

The  simple  iterative  procedure  in  which  the  separation  boundary  condition 
is  not  considered  failed  in  the  present  case  because,  whenever  the  free 
streamline  was  shifted,  its  length  changed  and  a  new  velocity  was  predicted 
in  order  to  satisfy  continuity.  The  resulting  potential  on  the  new  free 
surface  tended  to  vary  irratically  at  the  separation  point  leading  to  great 
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difficulty  in  generating  a  reasonable  flow  field  as  a  basis  for  Che  next 
trial-free  boundary.  However,  the  iteration  procedure  which  guarantees 
satisfaction  of  the  tangent  separation  condition  does  not  guarantee  that 
the  velocity  on  all  streamlines  will  be  the  same  because,  as  the  length  of 
the  streamline  changes,  the  average  velocity  on  the  streamline  also  changes, 
and  Chat  change  is  different  for  each  and  every  streamline.  If  the  flow  is 
close  to  axisymmetric ,  that  is,  if  the  cavity  shape  remains  essentially 
circular,  then  all  of  the  streamlines  tend  to  change  length  in  about  the 
same  amount  and  the  separation  tangency  condition  imposes  approximately  the 
same  constraints.  Accordingly,  the  solution  converges  very  nicely  in  that 
the  velocity  is  essentially  constant  over  the  entire  free  surface  (as  was 
found  with  the  circular-plane  and  some  of  the  elliptic-plate  pure-drag 
cases) .  In  the  next  section  the  detailed  results  of  these  calculations  are 
examined  and  hypotheses  are  developed  to  explain  the  behavior  and  potential 
solutions  for  the  problem.  Finally,  in  Section  6.3  examination  is  made  of 
alternative  formulations  which  might  also  solve  the  difficulty  encountered 


6.2  ACCURACY  OF  THE  SOLUTIONS 


6.2.1  Review  of  Tests  and  Relevant  Results 


As  indicated  in  the  introduction  to  Section  5,  the  separation  edge 
singularity  apparently  plays  a  crucial  role  In  the  accuracy  of  the  solutions 
for  various  three-dimensional  fully  cavitating  flows.  To  test  the  effect 
of  grid  refinement  near  the  separation  edge  four  additional  cases  were 
examined  for  the  circular  plate  pure  drag  flows  (cases  7-10  in  Table  5.1). 
The  grid  refinement  was  achieved  by  twice  halving  the  elements  adjacent  to 
the  separation  edge.  In  the  unrefined  cases  (4  and  5  in  Table  5.1)  the 
first  corner  node  from  separation  and  on  the  free  streamline  was  0.15  units 
downstream  of  the  circular  plate  (Ax^  *  0.15).  First-level  refinement  (see 
Figure  6.1)  used  Ax^  »  0.05;  second-level  refinement  used  Ax^  *  0.02. 

Both  long  (L  *  10)  and  short  (L  «  5)  cavities  were  studied.  Gross  changes 
in  a  ,  B/P  ,  C*  and  cavity  cross-section  were  less  than  one  percent  as  a 
result  of  refinement.  However,  significant  changes  were  observed  in  the 
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velocity  at  separation  and  the  consequent  variation  in  the  velocity  along 
a  given  streamline. 

Before  examining  these  results  a  brief  note  is  needed.  Larock  and 
Taylor  (1976)  point  out  that  with  the  finite  elements  used  in  their  paper 
(and  herein)  uniqueness  of  flow  velocities  is  obtained  by  averaging  all  the 
velocity  values  contributed  by  elements  which  contain  a  given  node.  As  an 
exception,  nodes  at  the  junction  of  a  solid  and  free  surface  can  be  con¬ 
sidered  to  have  in  common  only  the  solid  boundary  and  those  contributed  by 
elements  which  have  in  common  the  free  surface  boundary.  Larock  and  Taylor 
say  their  experience  suggests  that  one  should  average  only  the  velocities 
from  elements  having  a  free  surface  face.  Because  of  the  placement  of  mid¬ 
side  nodes  near  separation  at  the  quarterpoint  of  an  element's  side  (see 
Section  2. 2. 3. 3  or  Figure  (2.6),  the  x-  or  normal-components  of  the  velocity 
contributed  by  the  elements  with  solid  faces  and  contributed  by  elements 
with  free  surface  faces  are  exactly  equal.  However,  the  y-  and  z~  (or 
tangential-)  components  are  influenced  by  whether  or  not  the  element  has 
a  free  surface  or  solid  boundary.  The  present  tests  are  in  complete  agree¬ 
ment  with  the  Larock  and  Taylor  conjecture.  To  illustrate  these  points  the 
following  notation  is  used  in  subsequent  discussions:  the  velocity  at  separa¬ 
tion  contributed  by  elements  which  have  a  free  surface  face  is  denoted  by 

V_,  ;  the  velocity  contributed  by  elements  with  a  solid  face  is  denoted  by 
r 

Vg  ;  and  the  average  velocity  contributed  by  all  elements  at  a  node  is 

denoted  as  the  computed  velocity  V_  .  The  velocity  denoted  by  V  is  the 

free-stream  average  velocity  computed  by  summing  the  velocities  obtained  at 

the  cavity  centerline  on  each  streamline  and  dividing  by  the  number  of 

streamlines  (by  the  number  of  corner  nodes  in  the  y-direction;  see  App.  1). 

Figure  6.2  shows  the  effect  of  grid  refinement  near  separation  for  the 

circular-plate  pure-drag  cases.  Clearly,  the  accuracy  of  the  V_  represen- 

r 

tation  increases  dramatically  as  the  element  is  refined.  However,  V  ,  the 

w 

usual  average  velocity  at  a  node,  remains  a  constant  fraction  of  V  . 

r 

Therefore,  the  pressure  calculation  at  separation  will  remain  in  error  if  V 

w 

is  used  to  calculate  the  pressure  there.  However,  with  refined  elements 
near  separation  both  on  the  solid  plate  and  on  the  free  surface,  the  error 
in  the  pressure  or  drag  calculations  is  negligible.  If  the  finite  element 
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solution  were -exact  at  the  separation  point,  then  V.  would  equal  V  and 
both  would  equal  .  That  is  clearly  not  the  case.  Some  insight  can  be 
drawn  from  Figure  6.3  for  relatively  long  (L/P  ■  5)  and  short  (L/P  *  2.5) 
cavities  in  these  circular-plate  pure-drag  cases.  It  appears  that  as  W/P 
increases,  B/P  ,  i.e.,  the  width  of  the  cavity.  Increases  as  well.  Hence, 
the  curvature  of  a  streamline  near  separation  apparently  decreases.  One 
can  infer  that  the  "stress  at  separation"  should  decrease  as  W/P  increases. 
This  behavior  is  clearly  shown  in  Figure  6.3  where  both  Vg  and  Vc  increase 
relative  to  as  W/P  increases.  Therefore,  all  other  things  being  equal, 

the  separation  behavior  of  the  finite  element  method  becomes  less  critical 
as  the  wall  effects  decrease. 

Figures  6.4  and  6.5  illustrate  that  the  present  technique  is  a  very 

accurate  solution  along  each  individual  streamline  and,  further,  that  grid 

refinement  makes  a  significant  improvement  in  the  solution  behavior.  Figure 

6.4  shows  the  variation  of  the  free  stream  velocity  denoted  by  V__  with 

Fb 

distance  along  the  streamline.  This  variation  is  presented  as  a  percentage 
according  to  the  formula  AV/V(%)  ■  [(V_  -  V)/V]  x  100.  The  importance  of 

using  VF  at  separation  point  is  also  indicated.  For  example,  if  the 
averaged  velocity  at  the  separation  point  is  used  in  an  unrefined  calcula¬ 
tion,  the  apparent  velocity  error  at  separation  is  almost  22%.  However,  if 

V  is  used,  the  error  is  only  8%.  Similarly  the  maximum  error  away  from 
r 

separation  for  the  refined  circular  plate  case  is  less  than  1%.  For  the 
unrefined  case,  the  value  reaches  2-1/2%.  However,  for  most  of  the  stream¬ 
line  the  unrefined  and  refined  errors  are  essentially  indistinguishable.  A 
similar  plot  has  been  made  for  the  elliptic  plate  case  E2  and  is  given  in 
Figure  6.5.  This  figure  has  been  plotted  on  a  semi-log  scale  to  expand  the 
region  near  the  separation  point.  The  horizontal  axis  is  the  distance 
downstream  from  the  separation  point  plus  0.01  non-dimensional  units  so  that 
it  will  plot  on  a  log  scale.  It  is  important  to  notice  here  that  although 
no  specific  grid  refinement  was  undertaken  for  the  elliptic  case,  the  velocity 
error  is  less  than  2%  over  99%  of  the  streamline.  It  is  clear  then  that  the 
particular  method  used  here  is  a  highly  accurate  solution  along  each  stream¬ 
line. 

One  additional  experiment  was  also  performed  with  the  circular-plate 
pure-drag  case.  To  test  the  effect  of  increasing  the  number  of  elements  in 
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Che  transverse  direction  in  the  plane  of  section  A-A  (Figure  2.4)  ,  the 
number  of  elements  in  the  transverse  or  y-direction  (see  Appendix  1)  was 
doubled  for  one  pure-drag  circular-plate  case.  In  most  cases,  four  trans¬ 
verse  elements  were  used  in  the  quarter  flow  for  this  kind  of  computation. 

The  test  case  employed  8  elements.  No  gross  parameter  or  shape  changes 
occurred;  however,  the  accuracy  with  which  the  no  flux  boundary  conditions 
were  satisfied  on  the  y  ■  0  and  z  -  10  planes  of  the  flow  was  increased 
by  one  order  of  magnitude.  This  result,  together  with  the  refinement  tests, 
indicates  that  substantial  improvement  of  flow  detail  can  be  achieved  by 
adding  elements,  but  little  gross  behavior  change  appears  to  occur,  at  least 
in  these  pure-drag,  slightly  three-dimensional  cases.  The  information  gained 
by  this  test  of  increasing  number  of  elements  in  the  transverse  direction 
was  put  to  use  in  the  computation  for  the  elliptic  plate  pure  drag  cases. 

In  case  El  of  Table  5.2,  eight  elements  were  used  in  the  transverse  or 
y-direction.  However,  the  elements  nodal  points  were  more  or  less  uniformly 
distributed  in  the  transverse  direction  as  can  be  seen  from  Figure  5.4  This 
produced,  in  general  a  maximum  velocity  normal  to  the  vertical  boundary  on 
the  left  of  Figure  5.4  which  was  one-half  of  the  vertical  velocity  through 
the  upper  boundary  of  Figure  5.4.  In  addition,  the  computation  converged 
rather  raggedly.  In  both  cases,  E2  and  E3,  a  new  geometry  was  used  for  the 
elliptic  cases.  The  element  nodal  points  were  distributed  such  that  in 
each  element  the  included  angle  between  lines  drawn  from  the  element  corner 
nodes  to  the  center  of  curvature  of  the  free  surface  of  that  element  was 
constant  for  all  elements.  As  can  be  seen  from  Figure  5.4,  this  tends  to 
bunch  the  elements  near  the  sharp  end  of  the  ellipse  and  leave  rather  large 
spaces  on  the  flatter,  lower  end  of  the  ellipse.  Interestingly,  the  rearrange¬ 
ment  of  the  geometry  inverted  the  relationships  between  the  velocities  through 
the  vertical  and  horizontal  faces  of  the  flow,  the  maximum  velocity  through 
the  vertical  face  now  being  twice  the  maximum  velocity  through  the  horizontal 
face.  This  is  a  strong  indication  that  grid  refinement  in  the  streamwise 
direction  and  transverse  grid  refinement  and  geometry  are  important  in  fully 
three-dimensional  calculations . 

As  indicated  in  Section  3.3,  the  present  iteration  procedure  cannot 
guarantee  that  the  pressure  will  be  the  same  on  all  streamlines .  However, 
the  procedure  does  guarantee  that  the  tangent  separation  condition  is 
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satisfied  to  a  high  degree  of  accuracy.  Figure  6.6  is  the  variation  of 
free  stream  velocity  V?s  across  the  end  of  the  cavity  for  a  typical 
circular  plate  and  the  three  elliptical  plate  pure  drag  cases.  A  plot  is 
made  here  in  terms  of  node  points  in  the  transverse  direction.  Thus,  the 
actual  location  of  the  individual  streamlines  must  be  deduced  by  referring 
to  the  diagrams  for  the  various  flows.  In  any  case,  the  plot  moves  from 
left  to  right  in  the  same  manner  as  the  diagrams  for  the  various  flows. 

For  the  circular-plate  cases,  case  No.  5  which  is  plotted  here  contains 
the  maximum  error  of  any.  That  error  does  not  exceed  0.7  percent.  Thus, 
in  view  of  the  previous  results,  the  circular  plate  cases  can  be  considered 
to  be  exact  numerical  solutions  of  the  three-dimensional  fully  cavitating 
flow.  The  error  in  the  transverse  direction  for  the  free  stream  velocity 
was  not  significantly  influenced  by  the  streamwise  grid  refinement  that  was 
used  in  some  of  the  circular-plate  cases . 

The  elliptic  plate  pure  drag  cases  are  much  more  instructive  in  regards 
to  the  behavior  of  the  solution  because  the  streamlines  at  the  pointed  end 
of  the  ellipse  tend  to  be  much  more  highly  curved  than  those  at  the  flat  or 
middle  edge  of  the  ellipse.  The  response  of  the  computation  to  this  "stress 
at  separation"  condition  is  that  the  free  surface  velocity  is  much  lower 
at  the  highly  curved  edge  where  the  computation  is  unable  to  maintain  the 
tangent  separation  condition  at  a  higher  velocity.  Interestingly,  case  El 
which  uses  the  old  or  uniformly  distributed  geometry  for  the  elliptic  flow, 
behaves  relatively  more  uniformly  than  cases  E2  and  E3  where  the  method 
tends  to  break  down  on  the  smooth  edge  of  the  flow  where  there  is  very  little 
transverse  resolution.  This  is  consistent  with  our  earlier  observations 
about  the  flow  velocity  through  the  horizontal  and  vertical  no-flux  boundaries 
of  the  flow  in  the  elliptic  case. 

The  error  between  the  maximum  and  minimum  free  stream  velocities  on 
individual  streamlines  is  now  more  significant  in  the  elliptic  plate  cases 
ranging  up  to  a  maximum  value  of  about  7%  for  case  E2.  This  corresponds 
approximately  to  a  14%  error  in  the  constant  pressure  condition  if  one  makes 
the  computation  through  the  Bernoulli  equation  2.3. 

The  lifting  flow  with  a  circular  plate  set  at  an  angle  of  attack  of 
60°  in  a  rectangular  water  tunnel  has  more  extreme  conditions  as  far  as  the 
curvature  of  the  free  streamlines  than  either  of  the  pure  drag  cases  studied. 


From  Figure  5.6  it  can  be  seen  that  the  upper  or  leading  edge  streamline 
has  an  extremely  small  radius  of  curvature  near  separation,  while  the  lower 
or  trailing  edge  streamline  has  a  very  gentle  curve.  Thus,  the  free  stream 
velocity  computed  by  the  present  technique  is  much  higher  along  the  stream¬ 
line  from  the  trailing  edge  where  the  velocity  gradient  and  rate  of  curva¬ 
ture  is  smaller  so  that  the  tangent  separation  condition  can  be  more  easily 
satisfied.  If  streamwise  grid  refinement  were  the  essential  cure  for  the 
present  difficulties,  then  making  the  elements  near  the  leading  edge  separa¬ 
tion  should  decrease  the  relative  error  between  the  two  zones.  Thus,  one 
would  expect  to  more  equally  satisfy  the  tangent  separation  condition  and 
achieve  more  nearly  the  same  velocity  on  each  streamline.  To  pursue  this 
idea,  the  first  column  of  elements  of  separation  edge  were  made  smaller  at 
the  leading  edge  and  gradually  increased  as  one  moved  around  the  flow 
toward  the  trailing  edge  as  shown  in  Figure  6.7.  Though  this  refinement 
improved  the  convergence  characteristics  of  the  solution  and  produced  more 
rapidly  a  more  stable  solution,  no  substantial  improvement  was  observed  in 
the  error  which  was  approximately  10%  in  the  free  stream  velocity  on  the 
lower  versus  the  upper  streamline.  In  view  of  the  results  obtained  for  the 
cricular  plate  and  elliptic  plate  pure  drag  cases  and  the  result  obtained 
here,  it  would  appear  that  grid  refinement  should  be  made  in  the  transverse 
direction  and  that  this  would  lead  along  with  grid  refinement  in  the  stream- 
wise  direction  to  a  better  ability  for  the  computational  technique  to  meet 
both  tangent  separation  and  constant  velocity  (or  constant  pressure)  condi¬ 
tions  on  the  free  surface.  No  attempt  was  made  to  use  transverse  grid 
refinement  in  the  present  case  because  of  the  cost  of  the  computation. 

It  is  worth  noting  in  passing  that  during  the  course  of  the  iteration 
for  the  lifting  flow,  the  flow  rate  from  the  downstream  end  of  the  cavity 
generally  converges  to  its  final  value  which  is  very  close  to  the  inflow 
value  upstream  after  two  or  three  moves  of  the  free  surface.  The  flow  rate 
again  appears  not  to  be  a  good  criterion  for  measuring  convergence  because 
the  stream  surface  continues  to  move  significantly  until  the  fifth  or 
sixth  move.  Good  values  for  the  lift  and  drag  coefficient  are  also  obtained 
after  just  three  moves  indicating  that  for  gross  values  and  general  results 
is  not  necessary  to  iterate  to  a  final  solution. 
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One  final  test  was  made  with  a  lifting  flow.  The  tangent  separation 
condition  was  required  to  be  satisfied  only  at  the  leading  edge  of  the  cir¬ 
cular  plate.  The  separation  velocity  potential  was  then  adjusted  at  all 
other  streamline  in  an  attempt  to  create  a  free  surface  with  a  single 
constant  velocity.  The  solution  was  carried  through  six  free  surface  shifts. 
The  results  were  not  satisfactory  because  the  separation  velocity  at  the 
trailing  edge  exceeded  the  free  stream  velocity  and  the  cavity  surface 
showed  the  signs  of  drifting  experienced  by  Larock  and  Taylor  (1976) . 

6.2.2  Potential  Means  to  Reduce  Solution  Errors 

The  above  review  makes  it  clear  that  the  present  technique  is  satis¬ 
factory  only  in  cases  in  which  the  curvature  of  individual  streamlines  as 
they  separate  from  the  solid  body  is  essentially  the  same  over  the  entire 
separation  edge.  In  conditions  when  this  is  not  true  (for  example,  in  the 
lifting  flow  studied  here) ,  the  mixed  boundary-condition  singularity  at 
the  separation  point  causes  the  solution  to  be  inaccurate  and  to  make  it 
difficult  to  satisfy  both  the  constant  pressure  (or  constant  velocity)  con¬ 
dition  on  the  free  surface  and  the  tangent  separation  condition.  Also, 
mesh  refinement  both  in  the  streamwise  and  transverse  directions  tends  to 
reduce  the  computational  error.  There  appear  to  be  three  specific  things 
that  could  be  done  to  tackle  this  problem. 

The  first  and  most  direct  is  simply  to  increase  the  number  of  elements 
used  in  the  computation.  At  the  present  time  the  lifting  flow  utilizes  303 
elements  with  3,137  node  points.  Each  iteration  takes  5-1/2  minutes  of  CPU 
time  on  an  IBM  370/168  computer.  Since  approximately  thirty  iterations  are 
necessary  to  achieve  convergence  in  the  present  computation,  CPU  time  alone 
runs  in  excess  of  two  hours  per  computation.  On  the  other  hand,  Larock  and 
Taylor  (1976)  were  unable  to  examine  the  convergence  behavior  with  a  finer 
mesh  because  they  pointed  out  the  refinement  must  proceed  concurrently  in 
all  three  directions  or  else  some  elements  could  become  highly  distorted. 

In  fact,  in  the  present  case  the  elements  are  already  highly  distorted.  In 
particular,  because  of  the  relative  small  number  of  elements  in  the  trans¬ 
verse  direction,  the  elements  near  separation  could  be  viewed  as  long 
curved  pencils.  It  is  not  illogical  to  conclude  that  this  may  be  significant 
cause  of  the  difficulties  in  the  present  case.  Thus,  a  refined  mesh  seems 
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like  an  important  next  step  in  any  computations  for  lifting  flows  with  the 
present  method. 

In  order  to  accomplish  the  generation  of  a  refined  mesh  one  must 
examine  the  present  procedure  and  see  if  a  better  mesh  generator  can  be 
found.  In  the  present  mesh  generation  procedure  for  the  lifting  flow,  the 
plate  is  first  established  in  the  rectangular  tunnel  as  a  pure  drag  flow 
and  a  mesh  is  generated.  The  plate  is  then  tilted  to  the  desire  angle  of 
attack  and  the  elements  near  separation  are  carried  with  it.  The  elements 
at  separation,  particularly  around  the  leading  edge,  can  be  very  distorted 
in  the  streamwise  as  well  as  the  normal  direction,  although  they  are  not 
distorted  in  the  transverse  or  y  direction  particularly.  This  distor¬ 
tion  may  lead  to  difficulties  in  the  algebraic  equation  set  which  Is  gener¬ 
ated  to  solve  the  potential  flow  problem  or  may  produce  close  to  singular 
Jacobians  for  the  isoparametric  transformations. 

For  further  studies  and  as  a  second  step,  it  would  be  appropriate  to 
implement  a  mesh  generation  scheme  based,  for  example,  on  the  ideas  of 
Thompson,  et  al.  (1974) .  They  demonstrated  a  method  for  automatic  numerical 
generation  of  a  general  curvilinear  coordinate  system  with  coordinte  lines 
coincident  with  all  boundaries  of  a  general  multiconnected  region  contain¬ 
ing  any  number  of  arbitrarily  shaped  bodies.  They  stated  that  no  restric¬ 
tions  are  placed  on  the  shape  of  the  boundaries  and  the  method  is  not 
restricted  to  two  dimensions,  however,  it  is  not  known  if  the  method  has 
actually  been  implemented  for  a  three  dimensional  flow.  While  the  genera¬ 
tion  of  a  numerical  mesh  would  in  itself  be  expensive,  the  smooth  shape  and 
regular  coordinate  system  produced  by  such  an  automatic  generation  scheme 
for  curvilinear  coordinates  would  allow  the  user  to  obtain  a  very  accurate 
and  well  formed  finite  element  net  by  using  various  points  on  the  curvi¬ 
linear  coordinates,  translating  them  back  to  Cartesian  coordinates  and 
establishing  those  as  the  corner  and  mid  nodes  of  the  finite  elements.  With 
such  a  system  and  consistent  with  the  available  computation  funding  and 
raw  power,  one  could  accomplish  the  generation  of  an  adequate,  refined  and 
well-formed  finite  element  mesh.  Under  these  conditions  it  would  be  appro¬ 
priate  to  use  the  variant  on  the  present  iterative  technique  which  is 
applied  in  the  final  example  for  the  lifting  flow.  Namely,  control  the 
tangent  separation  condition  at  a  single  point  on  the  separation  edge  and 
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otherwise  adjust  the  separation  velocity  potential  to  produce  a  single 
constant  velocity  on  all  streamlines . 

A  third  and  final  alternative  would  be  to  follow  on  the  discussion 
in  Section  3.3  and  to  attempt  to  develop  a  formal  analytic  solution  for  the 
velocity  potential  in  the  zone  near  the  separation  edge.  This  would  have 
to  be  an  analytic  three-dimensional  solution  in  either  physical  or  iso¬ 
parametric  space.  If  it  could  be  found,  it  could  be  used  to  remove  singu¬ 
lar  behavior  from  the  finite  element  portion  of  the  problem  and  thereby 
allow  the  finite  element  solution  to  accurately  reproduce  the  regular  part 
of  the  flow.  The  development  of  such  an  analytical  patch  appears  to  be  a 
formidable  task  and  has  not  been  attempted  in  the  present  investigation 
(see,  however,  Strang  and  Fix,  1973,  for  an  outline  of  the  general  approach 
for  two-dimensions) . 

6.3  ALTERNATIVE  FORMULATIONS 

The  method  outlined  in  this  work  employs  three-dimensional,  20  node 
quadratic  finite  elements  and  a  trial-free-boundary  technique  to  solve  high 
speed  free  surface  flows.  As  noted  in  Section  4.2.1,  the  present  finite 
element  formulation  when  applied  to  a  two-dimensional  plane  flow  is  consider¬ 
ably  more  efficient  than  the  finite  difference  calculations  procedure  out¬ 
lined  by  Mogel  and  Street  (1974).  In  addition,  the  finite  element  formula¬ 
tion  meets  the  crucial  requirement  of  being  able  to  adjust  the  node  points 
of  the  element  grid  to  coincide  with  the  curved  and  irregular  boundaries  of 
the  free  surface.  It  has  been  made  clear  above,  however,  that  considerable 
grid  refinement  will  be  needed  to  make  the  present  method  an  efficient  and 
satisfactory  one  for  highly  three-dimensional  flows.  The  question  arises 
then  as  to  whether  or  not  there  is  a  better  alternative  formulation  which 
would  lead  to  a  successful  solution  of  nonlinear  fully-cavitating  flows. 

Two  alternative  formulations  come  easily  to  mind  and  are  discussed  below. 

The  first  is  the  use  of  continuous  finite  elements  with  the  hope  that 

these  elements  would  provide  higher  accuracy  at  the  separation  edge.  The 
second  method  is  the  boundary  integral  equation  method  in  which  it  is  neces¬ 
sary  to  solve  for  nodal  unknowns  only  on  the  surface  and  not  in  the  interior 
of  the  flow. 
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6.3.1  Other  Finite  Elements 

The  elements  used  In  the  present  computations  are  C°  elements,  that 
is,  they  insure  the  continuity  of  the  solution  function  4)  across  boundaries 
of  elements.  However,  the  first  derivatives  are  not  necessarily  continuous. 
The  present  elements  utilize  quadratic  interpolation  functions.  One  can 
define  C°  finite  elements  with  cubic  interpolation  functions.  However, 
there  is  no  inherent  gain  by  using  cubic  interpolation  functions  over  simply 
reducing  the  size  of  the  quadratic  interpolation  function  elements.  Accord¬ 
ingly,  it  does  not  seem  profitable  to  define  and  use  new  C°  elements. 

On  the  other  hand,  one  is  tempted  to  define  cubic  elements  using 

Hermite  polynomials  as  described,  for  example,  in  Pinder  and  Gray  (1977). 
Hermite  polynomials  with  an  isoparametric  approach  have  been  applied  in  two- 
dimensions.  However,  the  process  is  inherently  much  more  difficult  than 
that  used  here  and  it  is  not  clear  that  the  cost  of  using  the  Hermite  ele¬ 
ments  will  not  be  considerably  greater  than  using  the  present  C°  quadratic 
elements.  The  major  advantage  of  the  Hermitian  elements  is  that  the  unknowns 
now  involved  include  the  unknown  velocity  potential  <|>  and  its  first  deriva¬ 
tives  (i.e.,  the  flow  velocities)  at  the  comer  node-points  of  the  element. 
This  means,  for  example,  that  it  would  be  possible  to  specify  both  the 
velocity  potential  and  the  normal  derivative  of  the  velocity  potential  at 
the  separation  point  so  as  to  satisfy  both  the  constant  pressure  boundary 
condition  and  the  tangent  separation  condition  simultaneously.  Thus, 
Hermitian  elements  would  appear  to  be  attractive.  However,  Strang  and  Fix 
(1973)  have  conducted  a  numerical  experiment  on  the  behavior  of  various 
elements  in  the  neighborhood  of  a  mixed  boundary  condition  singularity  of 
the  type  encountered  at  separation  in  a  free  surface  flow.  Their  test 
indicate  that  the  best  accuracy  can  be  obtained  by  using  analytic  patches 
in  the  neighborhood  of  the  singularity  to  remove  the  singularity  so  that 
the  finite  element  does  not  have  to  cope  with  it.  On  the  other  hand,  a 
startling  result  is  that  the  standard  Hermite  cubics  come  out  worse  than  the 
simplest  of  linear  elements  with  regard  to  error  in  the  neighborhood  of  the 
singularity.  Strang  and  Fix  point  out  that  the  cubics  are  too  smooth  to 
cope  with  the  singularity.  It  would  appear  then  that  the  Hermitian  elements 
do  not  offer,  at  least  at  this  time,  an  attractive  alternative  to  the  solu¬ 
tion  method  used  so  far. 
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6.3.2  The  Boundary  Integral  Equation  Method  (BIEM) 

Brebbia  and  Wrobel  (1979)  and  Larock  (1977)  have  discussed  the  applica¬ 
tion  of  the  boundary  integral  equation  method  to  fluid  flows.  In  the  BIEM, 
all  computations  deal  directly  with  the  domain  boundaries,  in  this  case  the 
tunnel  walls,  solid  bodies  and  cavity  surfaces,  rather  than  requiring  a 
solution  throughout  the  entire  flow  domain.  If  the  BIEM  is  combined  with 
surface  finite  elements  then  one  can  utilize  finite  elements  in  much  the 
same  way  as  they  are  used  in  the  normal  finite  element  method.  That  is,  the 
body,  cavity,  and  tunnel  wall  surfaces  are  represented  by  a  set  of  finite 
elements,  e.g.,  quadratic  isoparametric  elements,  and  formulation  of  the  prob¬ 
lem  produces  a  series  of  nodal  unknowns  on  these  surfaces  only  which  may  be 
solved  in  the  standard  manner.  One  obvious  advantage  of  the  BIEM  with  finite 
elements  is  a  potential  substantial  reduction  in  the  number  of  nodal  unknowns. 

Larock  (1977)  focused  his  attention  on  jet  and  cavity  flows  and  out¬ 
lined  the  procedure  that  could  be  used  for  a  simple  two-dimensional  jet  flow. 
The  procedure  does  not  differ  in  the  large  from  the  one  used  here,  in  that  a 
trial-free-boundary  technique  is  used,  a  potential  flow  problem  is  solved 
and  the  resulting  velocity  field  is  used  to  move  the  free  surface  to  a  new 
trial  location.  The  only  difference  lies  in  the  method  used  solved  for  the 
velocity  potential  and  velocity  used.  BIEM  does  offer  potential  savings  in 
computer  time.  However,  the  problems  which  cause  difficulties  in  the 
present  computation  apparently  still  exist  in  the  boundary  integral  equation 
method . 

Larock  (1977)  indicates  that  his  formulation  of  the  BIEM  allows  the 
tangent  separation  boundary  condition  to  be  violated  and  it  was  only  by 
adjusting  the  velocity  on  the  free  streamline  that  he  was  able  to  drive  the 
normal  velocity  at  separation  to  zero.  This,  unfortunately,  is  precisely 
the  same  process  that  has  been  used  in  the  present  computation.  Therefore, 
it  would  appear  that  at  the  present  stage  the  boundary  integral  equation 
technique  has  the  same  failing  as  the  finite  element  technique  which  is 
implemented  above.  Accordingly,  it  is  probable  the  BIEM  can  be  used  only 
when  the  flow  field  is  essentially  axisymmetric  as  in  the  present  finite 
element  case  or  when  sufficient  grid  refinement  is  made  in  the  neighborhood 
of  separation.  Thus,  while  the  boundary  integral  equation  technique  should 
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probably  be  further  investigated  for  three-dimensional  fully  cavitating 
flow,  one  must  not  expect  it  at  the  moment  to  yeild  more  than  a  slightly 
more  efficient  computational  technique  and  the  difficulties  encountered 
in  the  present  investigation  will  not  be  instantly  overcome  by  shifting  to 
this  new  technique. 
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7.  SUMMARY  AND  FUTURE  WORK 


Based  on  potential  flow  theory,  a  finite  element  model  was  developed 
for  simulation  of  three-dimensional,  fully-cavitating  flows.  A  Riabouchinsky 
image  model  of  the  flow  was  used.  The  trial-free-boundary  approach  was  used 
and  effectively  reduced  the  nonlinear  potential  flow  problem  to  a  linear  one. 
However,  because  the  free  surface  shape,  as  well  as  the  velocity  potential, 
is  unknown,  iterative  solution  of  the  linear  problem  is  required. 

Larock  and  Taylor  (1976)  solved  a  three-dimensional  jet  flow  under 
the  Influence  of  gravity  by  use  of  a  method  and  computer  program  similar  to 
those  employed  herein.  However,  the  present  study  is  the  first  to  deal  with 
a  cavity  flow  and  the  first  to  conduct  a  detailed  examination  of  the  ability 
of  the  finite  element  method  to  accurately  represent  separated  flows.  This 
study  has  established  the  bases  for  evaluating  the  FEM  and  the  separation 
edge  problems.  The  explicit  role  of  the  free  surface  as  a  characteristic 
surface  of  a  quasi-linear  partial  differential  equation  was  demonstrated. 
While  this  result  is  of  no  interest  in  analytic  solution  approaches,  it  is 
crucial  for  numerical  approaches  and  is  at  the  core  of  the  separation  edge 
problems. 

A  unique  weighting  scheme  for  adjusting  the  free  surface  velocity 
potential  to  meet  the  required  tangent  separation  condition  was  devised. 

The  scheme  works  exceedingly  well  for  slightly  three-dimensional  flows. 
However,  in  highly  three-dimensional  flows  it  is  probable  that  both  a  very 
refined  grid  and  a  different  weighting  scheme  will  need  to  be  used. 

Both  pure  drag  and  lifting  flow  computations  were  made.  For  pure  drag, 
circular  and  elliptic  plates  in  a  water  tunnel  were  studied.  The  iterative 
procedure  was  stable  and  convergent.  The  solution  accuracy  was  excellent 
for  the  circular  plate  cases  and  acceptable  for  the  elliptic  plate  cases. 

One  lifting  flow  configuration  with  a  circular  plate  in  a  water  tunnel  was 
computed.  While  a  stable  and  converagent  solution  was  again  obtained,  the 
constant  pressure  condition  on  the  free  surface  was  not  accurately  satisfied. 

An  extensive  review  of  the  performance  of  the  solution  was  conducted. 
Several  conclusions  were  drawn.  First,  the  present  iterative  scheme  does 
not  guarantee  satisfaction  of  the  constant  pressure  condition  transversely 


across  Che  free  surface.  However,  this  is  of  concern  only  in  highly  three- 
dimensional  flows;  e.g.,  Che  pure  drag  solutions  studied  gave  quite  accu¬ 
rate  answers,  but  the  lifting  flow  solution  did  not.  Second,  grid  refine¬ 
ment  (reduction  of  element  size)  in  the  streamwlse  direction  and  near 
separation  markedly  improves  the  local  quality  of  the  solution,  but  has  little 
effect  on  global  results,  i.e.,  a  ,  C*  ,  etc.  Third,  grid  refinement  in  the 
transverse  direction  (perpendicular  to  the  free  surface  streamlines,  but  in 
the  free  surface)  will  be  necessary  to  achieve  accurate  results  in  highly 
three-dimensional  flows.  (A  relatively  coarse  grid  has  been  used  so  far 
because  of  cost.) 

Among  the  alternative  strategies  for  remedying  the  present  problems, 
only  grid  refinement  or  analytic  patches  to  remove  the  singular  nature  of 
in  the  neighborhood  of  the  separation  edge  offer  direct  solution  of  the 
difficulties.  While  the  boundary  Integral  equation  method  (BIEM)  may  reduce 
computation  cost  and  the  difficulties  associated  with  generating  a  "good" 
finite  element  net,  the  BIEM  appears  to  suffer  the  same  weaknesses  at 
separation  as  the  present  method. 

The  simplest  strategy  for  future  study  appears  to  be  to  generate  a 
highly  refined  FE  mesh  and  then  to  solve  again  the  present  lifting  flow  case. 
Use  of  the  BIEM  and  grid  refinement  is  the  next  most  difficult  strategy. 
Finally,  development  of  analytic  solutions  for  the  three-dimensional  flow  in 
the  neighborhood  of  separation  and  the  use  of  such  solutions  with  isopara¬ 
metric  elements  are  possible  in  principle,  but  are  well  beyond  the  present 
state-of-the-art . 
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APPENDIX  1 


MESH  GENERATION  AND  ELEMENT  ORGANIZATION 


One  of  Che  most  time  consuming  tasks  in  Che  FEM  is  setting  up  the 
geometric  input  for  a  not-so-regular  domain,  especially  in  three  dimensions. 
However,  by  use  of  a  semi-automatic  method,  the  geometric  input  for  this 
three  dimensional  cavitation  flow  (i.e.,  plates  in  a  rectangular  water  tunnel) 
can  be  set  up  relatively  painlessly.  The  flow  domain  is  divided  into  two 
parts,  namely,  a  rod  and  the  main  flow  section  (see  Figure  A. 1.1).  The 
tracking  of  the  node  numbers  is  accomplished  by  means  of  two  three-dimensional 
arrays,  namely,  MP  (ct,g,y)  for  the  main  flow  and  MPR  (a_ ,  8R>YR)  for  the 
rod.  Here,  a»  6  and  y  are  the  three  curvilinear  "axes"  for  the  main 
flow  as  shown  in  Figure  A. 1.1(a)  with  maximum  values  of  NROWS,  NCOLS  and 
NDEEP  respectively,  while  c^,  g^  and  y^  are  their  counterparts  for  the 
rod  section  as  shown  in  Figure  A. 1.1(b)  with  maximum  values  of  NROWR,  NCOLR 
and  NDEEPR  respectively.  Note  that  NDEEPR  is  function  of  the  value  of 
(counting  the  "0"  nodes).  Therefore,  with  a  given  set  of  values  of  a  , 
g  and  y  ,  one  can  pinpoint  the  node  number.  By  means  of  this  curvilinear 
axes  system,  the  free  surface  can  be  identified  atraight-forwardly ,  namely, 
a  *  1  and  g  >  NFREE  (Counting  from  separation  where  g  equals  NFREE) .  In 

the  present  configuration,  g  and  g  coincide  with  the  x-axis. 

K 


A.  1.1  Pure  Drag  Case 

Because  of  the  symmetry  In  the  flow  field,  only  \  of  the  flow  field 
need  be  considered  (refer  to  Figure  A. 1.2).  The  schematic  grid  is  first 
sketched  by  hand  such  that  elements  around  the  separation  are  smaller.  The 
size  of  the  rod  is  set  arbitrarily.  Then,  for  each  section  of  the  main  flow 
in  g-direction,  the  height  h  is  measured  together  with  h^  ,  h^  and  h^ 
(only  three  layers  are  shown  In  this  example) . 

The  ratios 


r 


1 


r 


2 


(A. 1.1) 
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and 


are  computed.  This  represents  the  distance  between  layers  in  the  a  direc¬ 
tion.  Using  ratios  r^  ,  r^  and  r^  ,  planes  for  each  y  (at  a  discrete 
value  of  9  )  are  discretized  accordingly.  The  sum  of  all  r  ratios  equals 

unity,  i.e. ,  £  r .  *  1. 

i=l  1 

A. 1.2  Lifting  Flow  Case 

Following  the  general  procedure  used  in  the  pure  drag  case,  one  first 
generates  a  pure-drag  half-flow  by  imaging  the  quarter  flow  as  shown  in 
Figure  A. 1.3.  Then  the  plate  is  tilted  to  the  desired  angle  of  attack  a  . 

The  grid  points  are  next  moved  proportionally  to  their  locations  with  respect 
to  the  separation  edge  and  the  fixed  tunnel  boundary  as  indicated  by  the 
dotted  lines  in  Figure  A. 1.3.  Figure  A. 1.4  shows  a  typical  arrangement  of 
the  grid  at  an  angle  of  attack  of  60  degrees. 

By  means  of  this  procedure,  one  can  generate  a  three  dimensional, 
lifting-flow  grid  relatively  painlessly.  However,  elements  around  the  lead¬ 
ing  edge  might  have  shapes  which  are  not  ''nice"  in  the  context  of  the  FEM. 

It  is,  therefore,  necessary  to  have  the  grid  plotted  and  examined.  If  the 
situation  arises  in  which  badly  distorted  elements  occur,  the  grid  is 
adjusted  semi -manually .  As  seen  in  Figure  A. 1.5,  it  is  desirable  to  have 
the  internal  angle  0^  smaller.  It  Is  done  by  pivoting  the  imaginary  line 
A-A  joining  points  T  and  B  at  some  point  and  rotating  the  line  to  position 
A'-A'.  Then  the  intermediate  nodes  are  moved  porportionally  on  the  line  A'-A'. 
The  position  of  the  pivot  governs  the  relative  movement  of  the  points  T  and  B. 
For  instance,  if  point  B  were  to  be  fixed,  the  pivot  is  set  at  that  point. 

The  input  sequences  and  the  program  listings  are  given  in  a  separate 
Program  Manual  (Ko  and  Street,  1979).  Therefore,  they  are  not  repeated  here. 
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APPENDIX  2 


LARGE  MATRIX  EQUATION  SOLVING  ALGORITHM 


As  mentioned  in  Section  2,  it  is  necessary  to  solve  a  large,  sparse 
system  of  algebraic  equations.  Due  to  the  fact  that  the  global  stiffness 
matrix,  in  Equation  2.25  is  positive  definite  and  symmetric,  one  has 

the  liberty  of  choosing  a  direct  method  over  an  iterative  method  or  vice 
versa.  The  choice  of  one  method  over  the  other,  generally,  depends  on  the 
availability  of  core  storage  and  one's  willingness  to  spend  computing 
time.  For  the  direct  method,  the  storage  requirement  depends  on  the  band¬ 
width  of  the  coefficient  matrix  of  the  system  and  on  a  clever  numbering 
scheme  of  elements  and  nodes  in  the  FEM  formulation.  In  the  present  study, 
a  direct  method  using  a  frontal  technique  with  Gaussian  Elimination  and  an 
iterative  technique  with  Conjugate  Gradient  method  were  tested.  The  itera¬ 
tive  Conjugate  Gradient  method  (ICGM)  as  implemented  here  enjoys  one  signifi¬ 
cant  advantage,  namely,  the  core  storage  is  only  a  fixed  multiple  of  the 
number  of  unknowns  and  is  independent  of  the  numbering  scheme.  This 
resulted  In  substantial  savings  in  core  storage  compared  to  the  direct 
method.  However,  the  overall  computing  time  is  longer  in  the  problem  tested 
here  for  the  ICGM. 

In  the  FEM,  the  global  stiffness  matrix  in  Equation  2.25  is 

usually  banded,  but  the  band  width  depends  on  the  type  of  elements  used  and 
the  numbering  scheme.  Often  it  is  this  global  matrix  together  with  fill-ins 
which  occur  during  application  of  direct  methods ,  which  require  more  core 
storage  than  is  available  economically  or,  often,  actually.  However,  K 


(e) 


ij 


In 


need  not  be  stored.  It  is  the  elemental  stiffness  matrices  K.  . 

ij 

Eqaution  2.22  which  must  be  saved  along  with  the  nodal  connection  matrix  v  . 
Both  the  direct  and  the  iterative  methods  operate  from  this  same  data  base. 


A. 2.1  The  Direct  Method  by  Gaussian  Elimination 

The  solution  of  the  FEM  formulation  equation  set  by  Gaussian  Elimina¬ 
tion  is  well  known  (Irons,  1970  and  Zienkiewicz,  1971).  The  method  described 
here  appears  in  a  program  supplied  by  C.  Taylor,  University  of  Swansea,  U.K. 
(see,  e.g.,  Larock  and  Taylor,  1976).  A  modified  version  of  this  program 
was  used  for  the  present  flow  problems. 
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The  elemental  stiffness  matrices  are  generated  and  stored  on  disk. 

They  are  read  into  core  one  at  a  time  filling  the  coefficient  matrix  in 
a  wedge  area  until  the  matrix  just  read  in  does  not  affect  the  first  row 
(or  equation)  in  the  "wedge"  (see  Fig.  A. 2.1),  which  moves  along  the  diagonal 
and  covers  the  upper  half  of  the  coefficients  present  in  the  global  stiffness 
matrix,  then  Gaussian  Elimination  is  performed  on  the  first  row.  The 
resulting  first  row  is  saved  on  another  disk  file  and  the  "wedge"  is  moved 
down  one  row.  More  elemental  matrices  are  then  read  in.  This  procedure  is 
repeated  n  times ,  where  n  is  the  total  number  of  nodes  in  the  problem. 
Next,  back  substitution  can  be  performed  by  recalling  the  rows  stored  on 
the  second  disk  file  one  at  a  time  in  inverse  order  of  storage  (last  in, 
first  out) .  The  core  storage  requirement  of  the  "wedge"  depends  on  the 
bandwidth  which  in  turn  depends  on  the  element  and  nodal  numbering  shcemes. 

No  pivoting  is  performed.  This  may  lead  to  loss  of  accuracy  in  some  problems, 
but  no  difficulties  were  encountered  in  the  present  cases  with  up  to  3137 
node  points. 

A. 2. 2  Iterative  Conjugate  Gradient  Method 

According  to  Concus,  et  al.  (1975),  the  Conjugate  Gradient  method 
was  proposed  by  Hestenes  and  Stiefel  in  1952.  It  was  not  until  the  early 
1970 's  that  it  was  found  to  be  highly  effective  as  an  iterative  procedure 
for  solving  large  sparse  systems  of  linear  equations.  The  method  is  based 
on  splitting  off  from  the  original  coefficient  matrix  a  symmetric,  positive- 
definite  matrix  that  is  more  easily  solved,  and  then  accelerating  the  iter¬ 
ation  using  conjugate  gradients.  The  convergence  of  the  CG  method  is  assured 
for  matrices  of  the  type  generated  by  the  FEM  formulation  described  above 
(cf.,  Concus,  et  al.  1975).  The  CG  method  is  described  briefly  below. 

Consider  the  system  of  equations 

Ax  ■  c  (A. 2.1) 

where  A  is  an  ixl,  symmetric,  positive-definite  matrix.  It  is  split 
such  that 

A  =■  M  -  N  (A. 2. 2) 
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where  M  is  positive-definite  and  symmetric  and  N  is  symmetric'.  Thus, 
from  (A. 2.1)  Mx  ■  Nx  +  c  .  It  is  assumed  that 


Mx  *  d 


can  be  solved  more  easily  than  (A. 2.1). 

Following  Concus,  et  al.  (1975),  the  CG  algorithm  can  be  stated  as 
follows: 

Let  x(0)  and  p^  ^  be  any  arbitrary  vectors.  For  k  =*  0,  1,  2,  ... 

(1)  Solve 


Mz 


(k) 


c  -  Ax 


(k) 


(2)  Compute 


2<k>  Mz(k> 
(k-l)T  Mz(k_1) 


,  k  >  1 


b  =  a 
o 


>(k)  -  z(k)  +btp(k-1) 

k 


(3)  Compute 


z<k>  Mz(k) 

(k)T  .  (k) 
p  Ap 


and 


(k+1)  (k)  .  (k) 

t  -  x  +  a^p 


The  above  algorithm  is  very  simple  and  very  easily  programmed.  The  maximum 
(k) 

value  of  a^P  in  step  3  above  is  used  as  the  maximum  error  between  iter¬ 
ations.  Instead  of  computing  the  right-hand  side  of  step  1  explicitly,  one 
could  use  the  following  recursive  form 


[c  -A*(lt+1)]  -  [=-««]-  Vp<k) 


(i.2.3) 


The  greatest  advantage  for  the  CG  method  is  that  one  has  complete 
liberty  in  choosing  the  matrix  M  in  (A. 2. 2).  If 


M  -  Diag  (A)  ,  (A. 2. 4) 

the  vector  z  can  be  very  easily  found.  The  matrix  A  is  used  only  once 
every  iteration  except  for  the  starting  iteration.  If  A  is  the  coefficient 
(global  stiffness)  matrix  K  in  (2.25) 


M  =  Diag  (K) 


(A.2.5) 


c  =  R  ,  and  the  final  x  is  the  solution  function  <P  in  (2.25).  In 
(k) 

computing  Ap  in  step  3  above,  one  need  only  use  the  elemental  stiffness 
matrices  K  from  (2.22)  as  follows.  For  each  element,  let  S  '  be 

defined  as  the  selection  matrix  of  dimension  n  x  l  with 


where 


.th  e  c(e) 

i  row  of  S 


=  yi(j) 


f  1  if  j  =  \ 

l  0  otherwise 


(A.2.6) 


(A.2.7) 


and  i(j)  are  all  distinct  for  i  »  1  ,  ...  ,  n  ;  then  it  can  be  easily 
seen  that 

Ap(k)  =  Kp(k)  =  Z  K(e)  S(e)  p(k)  (A. 2 

(e) 


However,  the  elemental  stiffness  matrices  need  to  be  preprocessed  to  incor¬ 
porate  the  Dirichlet  boundary  conditions  in  (2.18)  before  the  start  of  the 
CG  iterations. 
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Let  s  be  the  set  of  nodes  on  with  $  the  corresponding 
boundary  values.  Before  the  CG  method  can  be  applied  directly,  the 
elemental  matrices  (2.22)  have  to  be  preprocessed  as  follows: 


(i) 


(ii) 


(iii) 


(p) 

Set  ■  Kii  ^vi^  f°r  *  1  ,  ...»  n 

such  that  £  s 

For  i  =*  1,  2,  ...,  n,  such  that  v  e  s 

set  El  *  R  -  IC .  $(v.)  for  all  j  -  1,  ...»  n 
j  j  J 

such  that  (v^  e  s  and  ^  v^) 

Set  =0  for  all  i,  j  -  1,  2,  . . . ,  n 

such  that  i  ^  j  and  (v^  e  s  or  v  e  s) 


(A. 2. 9) 


(A. 2. 10) 


(A. 2. 11) 


Then  the  processed  elemental  and  the  global  matrix  equations  become 


K(e)  =  R(e) 


(A. 2. 12) 


K  <j>  -  K 


(A.  2. 13) 


respectively. 

It  can  be  easily  shown  that  the  submatrix  of  K  (with  the  known  values 
deleted)  is  essentially  a  principal  submatrix  of  K  and  is  therefore 
symmetric  and  positive-definite.  Accordingly,  the  CG  algorithm  can  be 
applied  to  solve  (A. 2. 13).  The  starting  vector  x  ,  of  course,  has  some 
known  values,  such  that 


x 


i 


for  all  i  £  S^ 


best  estimates 
or  simply,  0 


otherwise 


(A. 2. 14) 


72 


A. 2. 3  Numerical  Test 


The  CG  algorithm  and  the  direct  solution  method  were  tested  on  the 
IBM  370/168  of  the  Stanford  Linear  Accelerator  Center  Triplex  System.  For 
a  problem  with  1713  nodes,  304  elements  and  163  Dirichlet  boundary  nodes, 
the  estimated  bandwidth  of  K  was  160.  The  CG  algorithm  took  approximately 
0.5  sec  per  iteration.  Approximately  40  percent  of  this  time  was  used  in 
the  IBM  1/0  handling  routine,  30  percent  in  the  main  program,  and  30  percent 
in  the  system  routines  other  than  1/0.  In  the  case  of  double  precision 
arithmetic  (56  binary  bits) ,  the  solution  converged  to  machine  precision  in 
about  300  iterations.  Figure  A. 2. 2  shows  the  sequence  of  convergence. 

For  the  present  CG  algorithm,  one  must  save  five  vectors  of  the  size 
of  the  total  number  of  nodes  in  the  FEM  formulation.  In  this  test  case, 
the  CG  savings  in  core  storage  as  compared  to  the  direct  method  outlined 
in  Section  A. 2.1  was  40  percent.  However,  the  total  computing  time  for  CG 
as  compared  to  that  for  the  direct  method  was  approximately  1.5:1.  For 
another  larger  problem  with  2647  nodes,  496  elements,  303  Dirichlet  boundary 
nodes,  and  an  estimated  bandwidth  of  260  for  K  ,  the  CG  savings  in  core 
storage  is  approximately  60  percent. 
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Local  co-ordinates 
Mapped  Physical  Space 


Figure  2.3 
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Cartesian  co-ordinates 
Physical  Space 


The  isoparametric  finite  element 


Figure  2.4  Schematic  of  initial  cavity  flow  problem 


Element  at  separation  showing  quarter  point  location 


Definition  of  node  point  locations  for  periodic  cubic  spline 
interpolation 


Figure  3.2  Schematic  of  node  point  movement  for  lifting  flow 


Sequence  of  convergence  of  free  surface  for  Che  Cwo-diaensional 
cavity  flow 
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Figure  4.3 
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Figure  4. A  Demonstration  of  convergence  of  free-surface  shifting  algorithm 
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Figure  5.2  Variation  of  cavity  shape  with  tunnel  width  and  plate  diameter 
ratio.  P  •  2;  L  »  10;  cross-section  just  off  disk  and  at 
cavity  mid-section 


Figure  5.4  Variation  of  cavity  shape  with  tunnel  size  and  plate  diameter 
ratio  for  elliptic  plates.  Cross-section  just  off  disk  and 
at  cavity  mid-section 


Finite  element  grid  and  configuration  of  three-dimensional 
lifting  cavity  flow;  angle  of  attack  -  60° 


Figure  5.6  Final  cavity  shape  for  the  lifting  cavity  flow  with  circular 
plate 
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Figure  6.1 


Grid  refinement  at  separation 
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Variation  of  frees Cream  velocity  Vpg  (as  a  percentage 
difference  from  the  average  freestream  velocity)  for 
circular  plate  pure  drag  case 


Figure  6.5  Variation  of  freestream  velocity  VFS 


quarter  flow 


Figure  A. 1.4  Top  view  of  a  typical  grid  for  lifting  flow;  angle  of  attack 
of  60° 


Figure  A. 1.5  Schematic  of  semi-manual  grid  adjustment  at  the  separation 
for  the  lifting  flow 
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Figure  A. 2.1  Schematic  diagram  showing  the  data  structure  for  direct  method 
which  uses  Gaussian  Elimination 
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studied  were  fully  cavi tat ing  flow  past  flat  plates  in  a  water  tunnel.  Result 
are  given  for  pure  drag  flows  past  circular  and  elliptic  plates  and  a  lifting 
flow  past  a  circular  plate. 

Three-dimensional-,  20-node,  quadratic  isoparametric  finite  elements  are 
used.  The  locations  of  the  mid-side  nodes  on  the  free-surface  just  off  the  edf 
of  the  plate  are  found  to  have  significant  effects  on  the  separation 
condition  which  must  be  enforced  while  keeping  the  Jacobian  of  the  isopara¬ 
metric  transformation  nonsingular.  The  free  surface  is  shown  to  be  a  charac¬ 
teristic  surface  and  this  leads  to  the  development  of  a  weighting  scheme  for 
assigning  the  free  surface  potential  in  the  iterative  scheme  the  shifting  of 
the  free  surface.  A  new  three-dimensional  "local"  iterpolation  scheme  for  the 
location  of  the  mid-side  nodes  is  also  developed  for  the  lifting  flow  configu¬ 
ration. 

The  various  algorithms  were  tested  against  known  analytic  solutions  as 
well  as  published  numerical  and  experimental  results  in  two-dimensional  and 
axisymmetric  flows.  Results  from  several  pure-drag  three-dimensional  geom¬ 
etries  are  presented  and  the  cavity  shapes  and  the  tunnel  wall  effects  are 
examined.  Then  lifting  flow  results  are  given. 

^  Because  of  the  change  in  flow  boundary  conditions  at  the  separation  edge 
and  the  failure  of  the  FEM  to  resolve  these  conditions  accurately,  the 
ability  of  the  numerical  solution  to  maintain  a  constant  pressure  over  the 
entire  cavity  decreases  as  the  three  dimensionality  of  the  free  surface 
increases.  However,  the  present  procedure  produces  absolutely  stable  itera¬ 
tions  and  shows  no  sign  of  drifting  of  the  free  surface.  It  is  found  that 
satisfaction  of  a  tangent  separation  condition  of  the  free  surface  from  the 
flat  plate  body  is  crucial  for  the  stability  of  the  iterative  procedure. 

Grid  refinement  in  both  the  streamwise  and  transverse  directions  reduces  the 
computational  error.  While  free  surface  movement  between  iterations  is  a 
useful  convergence  criterion,  a  flow-rate  balance  between  upstream  and 
downstream  cross-sections  appears  not  to  be  a  good  criterion.  Finally, 
alternate  formulations  which  may  reduce  the  difficulties  encoun>«4:ed  are 
d i scussed . 
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